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ABSTRACT 


The  NPSAT1  satellite  uses  an  active  magnetic  torque  rod  system,  with  a 
magnetometer  for  attitude  determination,  to  maintain  3 -axis  stabilization,  with  a  slightly 
gravity  gradient  friendly  structure. 

This  thesis  will  examine  the  performance  of  three  combinations  of  programs  and 
simulation  models  for  the  NPSAT 1  satellite  attitude  control  system.  The  models  include 
a  magnetic  control  law  with  a  reduced  order  estimator  to  generate  torque  commands  to 
achieve  spacecraft  nadir  pointing  and  a  magnetic  rate  (Bdot)  control  law  to  reduce 
spacecraft  angular  rates.  The  perfonnances  of  two  Bdot  mode  switching  designs  are 
compared.  Also,  a  case  is  made  for  the  benefits  of  priming  the  system’s  reduced 
estimator  prior  to  mode  switching. 

All  of  the  control  methods  analyzed  appear  to  be  valid  control  methods  to  achieve 
three-axis  attitude  stabilization  using  only  magnetic  torquers  for  active  control.  The  most 
efficient  control  method  analyzed  incorporates  a  hand-off  method  from  a  magnetic  rate 
(Bdot)  control  loop  to  a  magnetic  control  loop.  The  results  of  this  analysis  indicates  that 
the  best  use  of  this  method  is  to  perform  the  Bdot  hand-off  following  the  achievement  of 
a  predetennined  combined  angular  rate. 
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EXECUTIVE  SUMMARY 


The  development  of  the  attitude  control  system  (ACS)  for  the  NPSAT1  satellite 
was  driven  primarily  by  one  object,  namely,  student  learning  through  analysis  and  design 
of  an  actual  spacecraft.  Another  concern  driving  the  NPSAT1  satellite’s  development 
effort  was  the  motivation  to  produce  a  small  inexpensive  research  satellite  composed  of 
COTS  (Commercial  Off  The  Shelf)  components  and  space  flight  experiments.  This 
naturally  led  toward  a  system  that  was  sensitive  to  the  anticipated  low  power  budget. 

During  the  effort  to  keep  the  above-mentioned  concerns  in  the  forefront  of  the 
attitude  control  system’s  design  process,  the  design  team  elected  a  control  system  that 
was  to  maintain  three-axis  stabilization  using  only  magnetic  torquers  and  magnetometers. 
Barry  Leonard  (Naval  Postgraduate  School  Adjunct  Professor)  produced  two  potential 
system  models.  Both  models  use  a  magnetic  control  law  and  a  reduced  estimator  to 
determine  the  spacecraft  states  and  achieve  three-axis  stabilized  nadir  pointing  using 
magnetic  torque  rods  and  a  magnetometer.  The  second  model  also  uses  an  added 
magnetic  rate  (Bdot)  control  law  as  suggested  by  the  Naval  Research  Laboratory  for 
initial  de-tumbling  of  the  spacecraft  after  deployment.  [Ref.  1] 

This  thesis  examines  the  performance  of  the  two  systems  as  stand-alone  models 
and  when  used  in  conjunction  with  the  Bdot  control  program.  One  application  of  the 
Bdot  control  loop  hands  off  attitude  control  after  a  predetennined  time  interval  and  a 
second  application  after  a  predetermined  combined  angular  rate  has  been  achieved.  After 
assigning  initial  conditions  and  establishing  a  test  strategy,  the  models  and  programs 
were  modified  to  include  data  capture,  storage  and  plotting  algorithms.  Due  to  the  need 
to  find  the  lowest  ACS  power  budget  possible,  in  order  to  allow  greater  power 
availability  for  research  payloads,  this  analysis  is  sensitive  to  the  trade-offs  between 
system  power  consumption  and  time  to  achieve  three-axis  attitude  control  after  launch 
and  deployment  or  after  a  loss  of  attitude  control  due  to  disturbance  forces.  The  ACS 
needs  to  achieve  and  maintain  stability  within  operational  requirements  (currently  nadir 
pointing  +/-  10  degrees),  yet  still  maintain  the  robustness  and  capability  to  recover  from  a 
loss  of  attitude.  The  results  from  each  of  three  combinations  of  control  techniques  is 


xv 


summarized  and  plotted.  Nine  sets  of  initial  conditions  were  tested  with  all  three  model 
and  program  combinations. 

All  of  the  control  methods  analyzed  appear  to  be  valid  control  methods  to  achieve 
three-axis  attitude  stabilization  using  only  magnetic  torquers  and  magnetometers.  The 
most  efficient  control  method  analyzed  appeared  to  be  the  model  that  incorporates  a 
hand-off  method  from  a  Bdot  control  loop  to  an  estimator  control  loop.  The  analysis 
results  indicate  that  the  best  use  of  this  method  is  to  perform  the  Bdot  hand-off  following 
the  achievement  of  a  predetennined  combined  angular  rate. 


I.  INTRODUCTION 


The  attitude  of  a  satellite  can  be  controlled  by  the  interaction  between  the  earth’s 
magnetic  field  and  a  magnetic  moment  produced  within  the  spacecraft.  Many  attitude 
control  systems  have  been  designed  using  magnetic  torquers  to  supplement  other  active 
or  passive  attitude  control  methods. 

The  purpose  of  the  NPSAT1  ACS  is  to  prove  a  satellite’s  attitude  can  be 
maintained  solely  with  magnetic  torquers  and  a  3-axis  magnetometer.  Magnetic  torquers 
rarely  fail  on  orbit  and  multitudes  of  satellites  have  magnetic  torque  rods  supplementing 
their  control  systems.  The  attitude  control  capability  of  the  magnetic  torquers,  however, 
has  limited  effectiveness  (+/-  several  degrees  of  pointing  accuracy).  One  possibility  for 
this  capability  is  a  low-power,  non-complex  safe  mode  for  a  more  complex  ACS.  This 
research  is  helpful  though  to  find  the  limits  of  technology  necessary  to  make  this  a  viable 
alternative  method.  The  possibility  of  a  low  cost  solution  to  attitude  control  that  could 
salvage  a  mission  after  complex  mechanical  failures  occur  would  be  a  substantial 
motivation  for  research  to  find  such  a  method. 

The  NPSAT1  satellite  will  use  an  active  magnetic  torque  rod  system  with  a  three- 
axis  magnetometer  for  attitude  detennination  to  maintain  3 -axis  stabilization  with  a 
slightly  gravity  gradient  friendly  structure.  This  thesis  will  examine  the  perfonnance  of 
three  combinations  of  programs  and  simulation  models  that  may  help  start  providing  the 
evidence  that  such  a  system  is  feasible  and  viable  for  the  NPSAT1  satellite.  Graphs  and 
tables  will  be  included  to  augment  the  presentation  of  the  analysis  data  in  this  thesis. 
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II.  BACKGROUND  AND  HISTORICAL  PERSPECTIVE 


NPSatl  is  the  second  satellite  in  the  Naval  Postgraduate  School’s  Small  Satellite 
Design  Program.  The  purpose  of  the  Small  Satellite  Design  Program  is  twofold.  First,  it 
provides  hands-on  education  for  officer  students  in  the  disciplines  of  satellite  design, 
system  integration,  space  systems  operations,  and  software  development.  Second,  it 
provides  students  with  dedicated,  on-orbit  assets  for  experimentation  and  testing. 

NPSatl  is  an  incremental  improvement  over  the  Postgraduate  School’s  first  small 
satellite,  PANSAT.  NPSatl  is  a  three-axis  stabilized  spacecraft  equipped  primarily  with 
Commercial-Off-The-Shelf  (COTS)  hardware.  Because  of  its  three-axis  stabilization,  it 
is  anticipated  that  the  spacecraft  could  support  a  basic  imaging  capability  with  a  COTS 
digital  camera. 

Several  spacecraft  components  have  already  been  acquired  or  built,  and  several 
experiments  or  payloads  have  already  been  envisioned  for  the  mission.  The  solar  cells 
and  basic  structure  of  the  spacecraft  were  inherited  from  a  cancelled  Navy  program  at  no 
cost  to  the  Naval  Postgraduate  School.  The  primary  spacecraft  controller,  a  modified 
COTS  80386  PC- 104  board,  was  already  designed  and  built.  The  Digital  Camera  and  its 
interface  board  had  been  purchased  and  tested  with  the  PC- 104  Central  Processor.  The 
Naval  Research  Laboratory  has  committed  to  flying  their  Coherent  Electromagnetic 
Radio  Tomography  (CERTO)  experiment  onboard.  Dr.  Alan  Ross,  Navy  TEN CAP  Chair 
in  the  Space  Systems  Academic  Group,  has  indicated  a  desire  to  test  his  Triple 
Redundant  Modular  Computer  experiment  on  board  NPSAT1  as  well. 

The  basic  orbit  profile  for  NPSAT1  was  to  be  a  circular  Low  Earth  Orbit  (LEO) 
between  400  and  800  kilometers,  inclined  at  between  40  and  80  degrees  with  a  mission 
life  of  18  to  24  months. 

The  NPSatl  Attitude  Determination  and  Control  System  (ADCS)  was  designed  to 
provide  three-axis  stabilization.  The  spacecraft  is  gravity  gradient  friendly  with  the 
possibility  of  an  optional  telescoping  boom  for  added  stability.  Active  attitude  control  is 
achieved  using  three  magnetic  torque  rods.  Attitude  detennination  software  receives 
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inputs  from  a  three-axis  magnetometer  to  determine  spacecraft  attitude  and  rates.  Figure 
2.1  shows  multiple  views  of  the  NPSAT1  spacecraft  and  the  planned  solar  array  layout. 
Three-axis  stabilization  provides  a  stable  spacecraft  to  support  the  imagery  and  CERTO 
Missions  on  NPSatl .  The  cone  is  at  the  top  of  the  spacecraft  and  the  nadir  pointing  face 
is  the  bottom. 


Figure  2.1.  NPS ATI 

The  NPSAT1  spacecraft  will  be  deployed  from  an  Air  Force  Delta  class 
expendable  launch  vehicle.  The  spacecraft  will  be  deployed  from  an  EELV  secondary 
payload  adapter  (ESPA)  ring.  Figure  2.2  shows  the  NPSAT1  spacecraft  as  it  will  be 
configured  on  the  ESPA  ring. 
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Figure  2.2.  NPSAT1  attached  to  EELV  payload  adapter  (ESP A) 
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III.  STATEMENT  OF  PURPOSE 


A.  OBJECTIVES  OF  RESEARCH 

Initially,  two  objectives  were  sought,  including  sensitivity  analysis  of  a  magnetic 
control  law  program,  and  after  a  suggestion  from  the  Naval  Research  Laboratory  to  try  a 
magnetic  rate  (Bdot)  control  law  [Ref.  1],  a  benefit  analysis  of  including  a  magnetic  rate 
control  law  for  initial  angular  rate  reduction  of  the  NPSAT 1  spacecraft.  In  the  course  of 
conducting  the  analysis,  two  secondary  objectives  were  identified.  Analysis  was  then  to 
include  the  difference  in  system  perfonnance  between  switching  from  Bdot  control  to 
nadir-pointing  control  after  a  time  interval  or  after  acquisition  of  a  specified  reduced 
angular  rate.  The  remaining  secondary  objective  was  to  explore  the  benefits  of  mode 
switching  before  or  after  priming  the  systems  reduced  order  estimator  with  historical 
data. 

The  analysis  was  conducted  to  detennine  the  best  configuration  of  the  existing 
proposed  attitude  control  systems  for  use  on  the  NPSAT  1  spacecraft.  Two  similar 
attitude  control  systems  were  examined  in  three  different  configurations.  The  first 
system  consisted  of  two  MATLAB  programs  and  one  SIMULINK  model  without  a  Bdot 
control  loop.  The  second  system  consisted  of  one  MATLAB  program  and  one 
SIMULINK  model  with  a  timed  hand-off  from  an  integrated  Bdot  control  loop. 
Additionally,  the  second  system’s  reduced  estimator  received  priming  data  prior  to  hand- 
off.  The  third  configuration  consisted  of  an  independent  Bdot  control  program  which 
executed  a  control  hand-off  to  the  first  SIMULINK  model  without  priming  the  reduced 
estimator,  after  a  predetermined,  reduced,  combined  angular  rate  was  achieved.  Table 
3.A.1  shows  the  programs  included  in  each  system  configuration. 

The  research  focused  on  three  objectives.  These  included:  initial  condition 
sensitivity  to  all  of  the  three  program  configurations,  the  difference  between  including  the 
Bdot  or  no  Bdot,  and  the  difference  between  a  timed  Bdot  hand-off  or  an  angular  rate 
based  hand-off. 
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SIMULINK  model 

MATLAB  system  programs 

Configuration  1 

LAB17.mdl 

Lab15MOIData.m 

"Labi 7  model" 

Lab16Data.m 

Configuration  2 

NPSAT  lACSINTRVW.mdl 

NPSAT  lACSData.m 

"Combo  model" 

Configuration  3 

Labi  7.mdl 

Lab15MOIData.m 

"Lab17  hand-off  model" 

Lab16Data.ni 

nps.m 

Table  3.A.I.  Configuration  models  and  programs 


B.  APPROACH 

Each  of  the  two  SIMULINK  models  was  designed  to  generate  control  torque 
commands  to  three  magnetic  torque  rods  for  active  attitude  control.  Using  initial 
spacecraft  ephemeris  data  and  a  spherical  hannonic  model  (8th  order),  the  systems 
estimate  what  the  Earth’s  magnetic  field  should  be  at  its  current  location.  The  system 
then  computes  the  cross  product  of  the  expected  magnetic  field  with  the  actual  local 
magnetic  field  measured  by  a  magnetometer.  A  cross  product  of  the  two  magnetic 
vectors  is  used  as  an  approximation  of  a  magnetic  control  law  to  produce  the  first  three 
system  states  (approximate  euler  angles).  A  reduced  estimator  is  then  used  to 
approximate  the  other  three  system  states  (angular  rates).  These  six  states  are  fed  to  an 
actuator  control  law  to  generate  torques. 

The  baseline  model  for  the  comparison  was  the  SIMULINK  model 
‘NPSATACSIntRvw.mdT  and  its  associated  m-file,  ‘NPSATlACSDATA.m’  which 
performed  a  timed  Bdot  control  hand-off  as  previously  described.  [Ref.  2]  This  model  is 
hereafter  referred  to  as  the  “combo”  model  since  the  Bdot  control  loop  is  integrated 
within  the  SIMULINK  model.  The  combo  model,  whose  block  diagram  is  shown  in 
Figure  3.B.1,  was  used  for  the  analysis  baseline  after  analysis  indicated  that  configuration 
exhibited  the  lowest  power  consumption  of  the  three  configurations.  The  three 
configurations  of  systems  were  evaluated  for  time  response  and  power  consumption. 
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Figure  3.B.I.  Block  Diagram  for  “Combo”  model 


The  programs  and  models  for  each  of  the  three  configurations  are  included  in 
Appendix  A. 
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IV.  PROGRAM  SETUP  AND  PROCEDURE 


A.  INITIAL  CONDITIONS  AND  ASSUMPTIONS 

In  order  for  the  attitude  control  system  for  the  NPSAT 1  spacecraft  to  be  viable,  it 
must  be  able  to  handle  the  conditions  immediately  after  deployment  from  the  launch 
vehicle,  and  handle  the  maintenance  of  the  spacecraft’s  operational  pointing 
requirements.  Initial  conditions  and  assumptions  were  established.  A  representative 
number  of  test  cases  were  run  to  establish  the  capability  of  the  system  to  respond  to  those 
conditions.  Each  configuration  was  compared  against  the  others  to  establish  the  best 
candidate  for  the  final  attitude  control  system.  The  primary  measures  of  effectiveness 
were  time  to  reach  steady  state  and  system  power  consumption. 

The  assumptions  and  initial  conditions  used  during  the  simulations  were  chosen  to 
approximate  as  closely  as  possible,  the  expected  final  characteristics  of  the  satellite  and 
its  anticipated  orbit.  The  minimum  allowable  altitude  due  to  aerodynamic  drag  and 
aerodynamic  disturbance  torques  for  the  spacecraft  is  400km.  The  expected  orbit  is 
circular  at  600km.  As  Figure  4.A.  1  shows,  the  three  30  AmA2  torque  rods  can  easily 
counter  solar,  gravity  gradient,  or  aerodynamic  torques  encountered  at  600km.  The 
generated  torques  are  greater  than  an  order  of  magnitude  more  powerful  than  the 
anticipated  disturbance  forces. 
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Torque  vs  Altitude 


—  Solar  Roll 

■  —  Solar  Pitch 
Solar  Yaw 

—  Grav  Grad  Roll 
— -Grav  Grad  Pitch 

Grav  Grad  Yaw 
— Aero  Roll 
Aero  Pitch 

—  Aero  Yaw 

—  30AmA2  Torq  rod  Min 
30AmA2  Torq  rod  Max 


Figure  4. A.  1 .  Torque  vs.  Altitude. 


Initial  satellite  angular  rates  were  based  on  the  expected  possible  tip-off  rates 
encountered  by  the  launch  vehicle  at  the  moment  of  deployment.  It  was  assumed  that  +/- 
0.1  rad/sec  was  a  conservative  range  of  possible  rates  for  spacecraft  pitch  and  roll.  Since 
the  yaw  rate  induced  by  the  launcher  to  the  spacecraft  is  not  expected  to  be  as  great  as  the 
pitch  and  roll  rates,  it  was  given  a  range  of  +/-  0.01  rad/sec.  For  simplification  and 
standardization  of  initial  conditions  between  tests,  the  initial  pitch  and  roll  rates  were 
assigned  either  a  “high  rate”  magnitude  (+  0.1  rad/sec  or  -0.1  rad/sec)  or  a  “low  rate” 
magnitude  (+0.01  rad/sec  or  -0.01  rad/sec).  The  initial  yaw  rate  was  assigned  a  positive 
or  negative  “low  rate”  only.  Considering  all  the  possible  combinations  of  the  “high”  and 
“low”  angular  rates  resulted  in  32  different,  testable,  initial  angular  rate  combinations. 

For  simplification  of  notation,  the  angular  rates  were  assigned  a  three  letter  combination 
of  either  capital  or  lower  case  N’s  or  P’s,  representing  the  magnitude  of  each  rate  as  well 
as  its  sign  (negative  or  positive).  Example:  nNp  represents  an  initial  pitch,  roll  and  yaw 
rate  of -0.01  rad/sec,  -0.1  rad/sec  and  +0.01  rad/sec  respectively. 

Power  consumption,  as  well  as  time  response  analysis  of  the  various  models  to 
each  of  the  selected  test  cases  required  the  capture  of  data  from  1 1  diferent  program 
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variables.  The  variables  included  the  output  torques  for  each  of  the  three  torque  rods,  a 
sum  of  the  absolute  values  of  the  three  rods’  output,  the  Euler  angles  <f),  0,  and  V| /,  the 
angular  rates  (f)dot,  0dot,  and  \|/dot,  and  time. 


B.  EXPERIMENTAL  PROCEDURE 

1.  Correlation  Between  Initial  Conditions  and  Test  Cases 

The  Bdot  control  law  program  ‘nps.m’  and  the  combo  SIMULINK  model  were 
each  evaluated  for  all  32  expected,  representative,  initial  angular  rate  combinations.  The 
programs’  time  responses  represented  the  time  each  program  required  to  reduce  the 
angular  rates  to  a  combined  rate  of  less  than  1.5  degrees  per  second.  The  resulting  time 
responses  seemingly  had  no  correlation  with  any  particular  aspect  of  the  initial 
conditions.  Table  4.B.  1 . 1  shows  the  lack  of  apparent  correlation  between  the  time 
responses  for  the  Bdot  control  law  program  ‘nps.m’  and  the  combo  model  and  the  initial 
conditions. 


■ 

Time  req'd 

Time  req'd 

Init  Rates 

Init  Rates 

Init  Rates 

Bdot  (sec) 

Init  Rates 

Bdot  (seel 

mm 

20000 

Pnp 

100000 

9200 

Pnp 

650 

£2" _ 

30000 

npn 

110000 

„PPD _ 

20 

ns 

550 

30000 

npP 

110000 

3500 

nPP 

640 

H3S11 

30000 

m _ 

110000 

IMSf 

6700 

m _ 

3700 

NPn 

30000 

110000 

NPn 

6900 

njisnn 

950 

nPn 

30000 

PNn 

110000 

nPn 

8900 

PNn 

14000 

UlJMau 

40000 

nnn 

120000 

17  7~^EL1 

11000 

nnn 

32 

ITS 

50000 

120000 

E9 

3500 

EESRt 

6700 

Pnn 

60000 

pNp 

120000 

Pnn 

1100 

pNp 

8700 

Nnn 

70000 

PPp 

120000 

Nnn 

3550 

PPp 

11000 

70000 

Pnp 

140000 

C7S 

3800 

m _ 

12 

pPn 

80000 

nPp 

150000 

pPn 

3500 

LUPUi 

8900 

pNn 

80000 

E2TSISII 

170000 

pNn 

8700 

ES 

9044 

NNn 

80000 

nNn 

170000 

NNn 

11000 

nNn 

9000 

nm 

90000 

PPn 

180000 

1 SH 

380 

PPn 

11200 

90000 

£ee _ 

270000 

840 

£EE _ 

3000 

Table  4.B.I.I.  Time  to  Achieve  System  Steady  State. 


Except  where  all  of  the  initial  rate  magnitudes  were  “low”  (0.01  rad/sec),  in 

which  case  response  times  tended  to  also  be  relatively  low,  the  data  did  not  support 

linking  higher  initial  angular  rates  with  longer  system  settling  times,  as  was  expected. 
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Neither  did  the  results  suggest  any  correlation  between  sign  order  of  the  initial  angular 
rates  with  a  higher  or  lower  system  settling  time.  Instead,  the  resulting  data  seemed  to 
suggest  that  the  ‘nps.m’  program  had  certain  initial  positive  or  negative  angular  rate 
combinations  that  the  program  responded  to  poorly  with  respect  to  system  settling  times. 
Many  combinations  whose  sign  order  matched,  yet  magnitudes  varied,  produced  similar 
system  steady  state  time  responses,  regardless  of  the  difference  in  individual  angular  rate 
magnitudes. 

The  data  did  not  favor  a  spread  of  test  cases  based  on  initial  angular  rate 
combinations  as  a  means  to  guarantee  a  representative  spread  of  system  performance. 
Instead,  using  the  data  collected  by  running  the  combo  model  simulation  (configuration 
2),  the  32  combinations  of  initial  angular  rates  were  sorted  in  order  of  their  system  steady 
state  settling  times.  The  combo  model  was  chosen  as  the  standard  of  perfonnance  to 
measure  all  other  programs  or  combinations  of  programs  against  to  detennine  the  optimal 
system  configuration  for  the  full  range  of  expected  initial  conditions.  Two  pieces  of  data 
of  specific  interest  in  the  analysis  was  the  total  system  torque  rod  output,  used  to 
calculate  an  energy  index,  and  the  time  to  reach  the  system  steady  state  for  each  test  case. 
Figure  4.B.  1 . 1  shows  an  early  plot  of  instantaneous  and  cumulative  torque  rod  power. 
Note  that  the  torque  rods  chosen  for  the  NPSAT1  spacecraft  are  limited  to  +/-  30  AmA2 
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Nine  test  cases  were  chosen  to  represent  a  spread  of  settling  times  demonstrated 
by  the  combo  model,  from  the  lowest  demonstrated  settling  time  to  the  highest.  The  nine 
sets  of  initial  angular  rates  included  “ppn,”  “NPn,”  “NNp,”  “pPn,”  “Pnp,”  “nnn,”  “pnp,” 
“nNp,”  and  “Ppp.” 

2.  Run  Program,  Capture  Data 

There  were  three  configurations  of  programs  and  models  that  were  run  for  each  of 
the  nine  test  cases.  For  each  iteration,  the  data  was  recorded  and  analyzed  to  determine 
the  power  required  by  the  system  to  reach  the  desired  system  steady  state.  The  power 
required  is  represented  in  the  form  of  an  energy  index.  The  torque  rod  output  is 
measured  in  AmA2.  If  the  output  is  integrated  over  the  length  of  time  to  reach  system 
steady  state,  the  result  is  an  “energy  index”  (AmA2  *  sec)  that  was  used  for  comparison 
against  the  other  three  programs  and  methods  that  were  run  with  the  same  initial 
conditions. 

In  order  to  capture  the  data  variables  of  interest  from  the  simulations  and  program 
runs,  the  programs  and  simulation  models  were  modified  to  include  output  blocks  to 
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generate  MATLAB  ‘.mat’  files  which  included  the  desired  data,  or  the  programs  were 
modified  to  generate  new  variables  which  could  be  saved  from  the  MATLAB  work 
space.  The  programs  ‘nps.m’  and  the  combo  model  were  each  run  for  all  32 
combinations  of  initial  angular  rates. 

The  Bdot  control  law  program  ‘nps.m’  required  modifications  of  the  MATLAB 
source  code  to  record  the  1 1  variables  from  each  program  run.  Additional  data 
manipulation  and  plotting  programs  were  created  for  further  data  analysis  and  to  plot  the 
analysis  results.  The  program  ‘grabdataset.m’  was  written  for  execution  following  every 
different  iteration  of  the  ‘nps.m’  program  to  clear  the  MATLAB  workspace  after  copying 
all  desired  variables  into  an  1 1  element  matrix.  Each  of  the  variables  was  sampled  at  a 
pre-assigned  change  in  time  or  dt  (sec)  until  the  simulation  reached  its  chosen  final 
combined  rates  goal.  The  value  chosen  for  all  of  the  simulations  was  dt  =  2.0  seconds. 
The  plotting  programs  ‘plotresult2.m’  and  ‘rodtotal  vstimc.m’  were  used  to  provide 
graphical  representation  of  the  data.  The  source  code  for  the  programs  ‘grabdataset.m,’ 
‘plotresult2.m’  and  ‘rodtotalvstime.m’  are  provided  in  Appendix  A. 

Both  SIMULINK  models  were  modified  to  capture  the  same  1 1  variables  that 
were  recorded  from  the  ‘nps.m’  program  runs.  Each  SIMULINK  model  was  run  with  the 
same  initial  conditions  used  in  the  iterations  that  were  run  with  the  Bdot  control  law 
program.  The  output  blocks  added  to  the  simulation  models  are  shown  in  Figure  4.B.2. 1 . 

3  rod  torques  Euier  ang|es  Angular  rates 


Figure  4.B.2. 1 .  Simulink  Output  Data  Capture  Blocks. 


16 


The  output  ‘.mat’  files  were  imported  into  the  MATLAB  work  space  for  data 
analysis  and  plotting  using  the  programs  ‘grabdataset.m,’  ‘create_variables.m,’ 
‘rodtotalvstime.m’  and  ‘plotresult2.m,’  which  are  provided  in  Appendix  A. 

The  program  ‘grabdataset.m’  was  written  for  manipulating  data  captured  from  a 
SIMULINK  simulation  run.  The  program  measured  and  plotted  the  instantaneous  and 
cumulative  total  of  commanded  rod  torques  and  the  time  required  to  achieve  a  3 -axis 
steady  state  pointing  accuracy  of  less  than  10  degrees.  This  program  also  solved  and 
plotted  the  area  under  the  instantaneous  rod  torque  total  curve,  producing  an  energy  index 
for  comparison  with  other  runs  using  different  initial  conditions. 

Eleven  variables  were  recorded  in  a  matrix  as  a  ‘.mat  file’  from  a  simulation  using 
one  of  two  modified  SIMULINK  simulations.  The  program  picked  off  the  following 
desired  variables:  the  torque  rod  outputs  in  AmA2  as  ‘rodl’,  ‘rod2’,  and  ‘rod3.’  The  total 
commanded  torque  from  the  three  rods  combined  was  ‘rodtot.’  The  angular  rates  were 
recorded  as  ‘Phd,’  ‘Thd’  and  ‘Psd,’  and  the  Euler  angles  were  recorded  as  ‘Ph,’  ‘Th’  and 
‘Ps.’ 

The  Bdot  MATLAB  program  output  data  and  SIMULINK  system  models’  output 
data  needed  to  be  merged  into  a  common  format  for  data  analysis.  The  program  ‘create 
variables. m’  was  written  for  capturing  data  from  the  NRL  program  ‘nps.m.’  The 
program  captured  1 1  variables  for  analysis  and  plotting  purposes.  Many  of  the  data 
parameters  were  manually  entered  into  the  program  code  from  the  MATLAB  workspace. 
Future  analysis  will  benefit  from  the  development  of  an  automated  data  capture  program 
that  combines  all  of  the  steps  used  to  manipulate  the  data  for  this  analysis.  The  resulting 
data  was  saved  from  the  MATLAB  workspace  to  a  data  folder.  The  ‘.mat’  file’s  name 
reflected  the  initial  conditions  used  to  produce  the  data. 

EXAMPLE  FILE  NAME:  ‘bdot_pNp_casel.mat.’  In  this  example,  pNp  stands 
for  ‘Phd’  =.01  rad/s,  ‘Thd’  =-.l  rad/s,  ‘Psd’  =.01  rad/s. 

The  program  ‘RodtotalvsTime.m’  was  written  to  plot  simulated  torque  rod  data 
for  analysis  of  simulations  supporting  NPSSAT1.  This  program  captured  a  row  of  data 
from  a  saved  1 1  row  matrix.  The  matrix  was  produced  from  one  of  two  modified 
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SIMULINK  simulations.  The  program  grabbed  torque  rod  totals  and  plotted  them  vs. 
time. 

3.  Run  Program  Handoff  Method,  Capture  Data 

The  program  ‘NPSSATlTESTCASE.m’  was  written  for  manipulating  data 
captured  from  MATLAB  programs  and  SIMULINK  runs  that  were  run  in  the  hand-off 
mode.  Since  each  program  ran  independently,  their  data  had  to  be  combined  for  analysis. 

This  program  computed  and  plotted  the  results  of  the  two  different  simulations’ 
data.  The  Bdot  control  law  program  ‘nps.m’  was  used  to  arrest  angular  rotation  rates  to 
the  predetermined  rate  goal  (less  than  or  equal  to  1.5  degrees/sec).  Then,  the  final  states 
of  that  simulation  were  used  for  the  initial  conditions  of  the  second  program,  which 
further  arrested  the  rotation  rates  and  achieved  a  predetermined  three-axis  nadir  pointing 
accuracy.  Data  from  the  program  ‘nps.m,’  was  used  to  provide  the  initial  conditions  for 
each  of  two  modified  SIMULINK  simulations.  The  data  stored  as  a  result  of  running  one 
of  the  two  follow-on  simulations  were  loaded  into  the  MATLAB  work  space,  with  the 
stored  data  from  the  Bdot  program  run.  The  ‘NPSSATlTESTCASE.m’  program 
essentially  combined  two  other  programs,  ‘grabdataset.m’  and  ‘create  variables.m.’ 

After  loading  the  desired  ‘.mat’  files  containing  data  matrices  into  the  MATLAB 
workspace,  the  matrix  sizes  were  used  to  edit  the  program’s  array  sizes.  The  data  are 
saved  as  a  single  ‘.mat’  file. 

The  ‘NPSSATlTESTCASE.m’  program  captured  torque  rod  total  and  time 
variables  for  analysis  and  plotting  purposes.  The  program  then  measured  and  plotted  the 
instantaneous  and  cumulative  total  of  commanded  rod  torques  and  the  time  required  to 
achieve  a  3-axis  steady  state  pointing  accuracy  of  less  than  10  degrees.  This  program 
computed  and  plotted  the  area  under  the  instantaneous  rod  total  curve  to  produce  an 
energy  index  for  comparison  with  other  runs  using  different  initial  conditions. 

4.  Run  Plotting  Programs 

Once  the  four  programs  or  combinations  of  programs  were  run  for  each  of  the 
nine  test  cases,  the  ability  to  use  and  manipulate  the  data  with  calculations  was  relatively 
straight  forward.  In  order  to  manipulate  the  data  from  both  the  independent  Bdot 
program  and  its  second  handoff  program,  the  data  had  to  be  merged  together  at  the  point 
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of  handoff,  which  occurred  at  the  instant  the  Bdot  program  reached  its  assigned  combined 
angular  rate  threshold.  The  resultant  data  was  then  used  as  if  it  were  merely  the  results  of 
one  of  the  individual  programs. 


Figure  4.B.4.I.  Angular  Rates  vs.  Time  (pPn). 


An  example  of  the  plotting  program  ‘plotresult2.m’  is  shown  in  Figure  4.B.4.I. 
The  initial  angular  rates  were  <f)dot  =  +0.01  rad/sec,  Odot  =  +0.1  rad/sec,  and  \|/dot  =  -0.01 
rad/sec  (pPn).  The  program  ‘nps.m’  was  assigned  a  combined  angular  rate  goal  of  1.5 
deg/sec.  The  program  arrested  the  angular  rates  until  the  sum  of  the  absolute  values  of 
the  three  angular  rates  dropped  below  the  assigned  threshold. 

Figure  4.B.4.2  shows  the  results  of  configuration  3  (on  the  left  side  of  the  figure), 
the  Bdot  control  handoff  method  to  the  ‘Labl7.mdl’  for  the  test  cases  “pPn”  (4>dot  = 
+0.01  rad/sec,  Odot  =  +0.1  rad/sec,  and  \|/dot  =  -0.01  rad/sec)  and  “NNp.”  (4>dot  =  -0.1 


19 


rad/sec,  0dot  =  -0. 1  rad/sec,  and  \|/dot  =  +0.01  rad/sec).  In  this  method,  the  SIMULINK 
model  ‘Labl7’  has  not  been  activated  until  the  Bdot  program  ‘nps.m’  has  finished 
arresting  the  angular  rates. 
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Figure  4.B.4.2.  Bdot  ‘nps.m’  handoff  to  ‘Labl7.mdl’. 


The  point  at  which  the  handoff  occurred  is  quite  distinguishable  in  these  cases, 
occurring  at  5,000  seconds  and  1,500  seconds  respectively.  The  graphical  representation 
of  the  data  in  this  case  was  instrumental  in  recognizing  a  difference  in  the  smoothness  of 
transitions  between  different  handoff  models.  Using  the  same  initial  conditions  for 
configuration  2,  the  internal  system  Bdot  control  handoff  to  the  estimator  in  the  combo 
SIMULINK  model  shows  a  much  smoother  transition  at  the  handoff  point.  As  Figure 
4.B.4.2  indicates,  the  sudden  spike  in  required  energy  from  the  torque  rods  does  not 
occur.  This  observation  proved  to  be  true  for  many  of  the  handoff  cases. 


20 


While  attempting  to  explain  why  there  was  a  difference  between  the 
configurations,  several  insights  were  gained  into  the  operations  of  the  models.  The 
combo  model,  essentially,  is  a  Bdot  loop  operating  for  a  period  of  time  within  the  model, 
then  handing  off  the  control  to  the  reduced  order  estimator.  The  Kalman  filter  within  the 
seemingly  idle  reduced  order  estimator  loop  of  the  model  was  in  fact  gathering  magnetic 
field  crossed  with  magnetometer  infonnation,  while  the  Bdot  loop  was  arresting  the 
initial  angular  rates.  Therefore,  when  the  control  was  handed  off  by  the  Bdot  loop  to  the 
reduced  order  estimator,  the  filter  already  had  better  knowledge  than  if  the  system  started 
gathering  magnetic  field  data  at  the  handoff  point.  The  right  side  of  Figure  4.B.4.2  shows 
the  handoff  which  includes  the  “better  knowledge.” 
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V.  RESULTS 


A.  OVERVIEW 

The  progression  of  analysis  objectives  began  with  sensitivity  analysis  of  the 
existing  ACS  model,  which  had  no  Bdot  control  law.  While  that  analysis  was  being 
conducted,  Naval  Research  Laboratory  personnel  recommended  using  a  magnetic  rate 
(Bdot)  control  law  to  reduce  the  satellite  angular  rates  to  a  lower  “manageable”  level 
prior  to  using  the  reduced  estimator  models.  [Ref.  1]  The  objective  of  showing  the 
difference  in  perfonnance  between  systems  with  or  without  the  Bdot  control  law  was 
started.  It  was  expected  that  the  systems  would  acquire  steady  state  nadir  pointing  faster 
if  the  initial  angular  rates  were  reduced  faster.  The  Bdot  program  ‘nps.m’  was  analyzed 
for  32  different  sets  of  initial  conditions  and  reduced  the  rates  much  faster  than  the 
reduced  estimator  control  system  without  Bdot.  The  performance  measurements  were 
based  on  the  time  required  for  the  system  to  achieve  nadir  poining  accuracies  of  +/-  10 
degrees. 

After  several  evolutionary  steps  of  the  ACS  system,  the  combo  model  was 
produced,  which  incorporated  the  Bdot  control  law  within  the  system  model.  The  power 
consumption  results  out-perfonned  all  previous  system  versions.  The  decision  to  use  the 
combo  model  as  a  benchmark  for  perfonnance  was  made.  The  combo  model, 
configuration  2  was  analyzed  for  all  32  initial  test  case  combinations  for  time  response. 
Curiously,  the  time  responses  for  the  combo  model  compared  to  the  other  configurations 
did  not  have  the  best  performance.  It  was  discovered  that  a  difference  in  time  response 
was  achieved  if  the  Bdot  hand-off  was  conducted  after  a  predetermined,  reduced  angular 
rate  was  achieved  instead  of  after  a  predetermined  time  interval.  The  decision  to  conduct 
analysis  in  support  of  three  objectives  was  then  made.  The  three  objectives  included: 
initial  condition  sensitivity  to  all  of  three  program  configurations,  the  difference  between 
Bdot  and  no  Bdot,  and  the  difference  between  timed  Bdot  hand-off  and  angular  rate 
based  Bdot  hand-off. 
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B. 


BDOT  RESULTS 


The  results  of  the  Bdot  program  runs  offered  no  insight  into  the  correlation 
between  initial  conditions  and  system  time  response.  Sorting  the  data  by  time  required  to 
reach  the  combined  angular  rate  goal  of  <=1.5  degrees/sec,  had  no  correlation  with  a 
similar  sorting  of  any  of  the  other  systems’  results.  There  is  no  doubt  that  the  Bdot 
control  law  program,  ‘nps.m,’  by  itself,  arrested  the  angular  rates  faster  than  any  other 
program  or  model.  (See  Figure  5.B.1) 
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Figure  5.B.I.  Bdot  Time  to  Reduce  Rates  to  <=  1.5  deg/sec  Combined. 


Unfortunately  the  program  does  not  offer  any  three-axis  pointing  capability. 
There  was,  however,  substantial  time-savings,  often  by  a  factor  of  ten  or  more  over  the 
ability  of  the  other  programs  to  achieve  partial  acquisition  by  angular  rate  reduction. 
This  performance  led  to  the  decision  to  analyze  the  affect  of  using  the  Bdot  program 
independently  until  a  specified  combined  angular  rate  goal  was  achieved.  The  final 
conditions,  including  angular  rates  and  Euler  position  angles  were  handed  off  to  the 
original  ‘Labl7.mdl’  as  initial  conditions,  (configuration  3). 

The  Bdot  control  law  more  efficiently  arrested  the  angular  rates,  compared  to 


configuration  1,  the  system  without  the  magnetic  control  law  (‘Labl7.mdl’).  The  Bdot 
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program  also  appeared  to  help  the  other  systems  most  efficiently  when  the  Bdot  program 
was  run  to  a  specific  rate  goal,  instead  of  for  a  specified  time.  Once  the  control  was 
handed  off  to  the  other  programs  or  models,  the  energy  levels  and  commanded  torques 
required  to  resume  the  control  effort  were  minimal. 


C.  CONFIGURATION  1  RESULTS 

After  conducting  the  nine  test  cases,  configuration  1  seemingly  had  the  poorest 
performance  of  the  three  configurations.  Compared  to  the  other  two  configurations,  to 
achieve  pointing  accuracies  of  +/-  10  degrees,  configuration  1  appeared  to  only  perfonn 
well  in  the  middle  ground  between  the  extremes  of  initial  conditions  (using  the  combo 
model  as  a  bench  mark).  The  best  cases  and  worst  cases  for  configuration  2  (the 
extremes),  however,  nearly  always  corresponded  to  poor  performance  from  this 
configuration. 

The  energy  indexes  for  this  configuration  tended  to  exceed  the  worst  performance 
of  the  other  two  methods.  The  gains  selected  for  this  model,  however,  have  been 
designed  to  be  adequate  for  the  entire  range  of  initial  conditions.  It  is  quite  likely  that  if  a 
Bdot  control  law  was  added  to  this  model,  and  the  Bdot  program  maintained  its  control 
until  a  minimal  residual  combined  angular  rate  existed,  then  the  gains  for  the 
‘Labl7.mdT  could  be  optimized  to  accommodate  the  low  initial  rate  conditions.  Table 
5.C.  1  shows  the  nine  case  studies  conducted.  The  order  of  the  initial  conditions  is  in 
increasing  system  time  order  for  the  reference  combo  model.  Figures  5.C.  1  through 
5.C.18  show  the  instantaneous  and  cumulative  power  consumption  results  for 
configuration  1. 
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Lab  17 

Energy  Index  [(Am*2nsec]x1e4] 

Init  Rates 

Lab  17 

200000 

119278 

MPn 

67000 

169540 

NNp 

92000 

5585031 

pPn 

60000 

66757 

Pnp 

90000 

2474116 

nnn 

30000 

1 

170000 

98606 

200000 

2677909 

IfEE _ 

40000 

153848 

Table  5.C.I.  Configuration  1  (‘Labl7.mdT)  test  results. 
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Figure  5.C.I.  “ppn”  Configuration  1:  Cumulative  Rod  Total  vs.  Time 
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Figure  5.C.2.  “ppn”  Configuration  1:  Instantaneous  Rod  Total  vs.  Time 
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Figure  5.C.3.  “NPn”  Configuration  1:  Cumulative  Rod  Total  vs.  Time 


Figure  5.C.4.  “NPn”  Configuration  1:  Instantaneous  Rod  Total  vs.  Time 


28 


NNp(Lab17)  area  Rod  Total  558503.0668  in  100000  sec 


Figure  5.C.5.  “NNp”  Configuration  1:  Cumulative  Rod  Total  vs.  Time 

NNp  (Lab17)  Rod  Total  vs.  Time 


Figure  5.C.6.  “NNp”  Configuration  1:  Instantaneous  Rod  Total  vs.  Time 
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pPn(Lab17)  area  Rod  Total  66756.538  in  70000  sec 


Time  (sec) 


Figure  5.C.7.  “pPn”  Configuration  1:  Cumulative  Rod  Total  vs.  Time 


pPn  (Lab17)  Rod  Total  vs.  Time 


Figure  5.C.8.  “pPn”  Configuration  1:  Instantaneous  Rod  Total  vs.  Time 


30 


Pnp(Lab17)  area  Rod  Total  247411.6284  ir  90000  sec 


Figure  5.C.9.  “Pnp”  Configuration  1:  Cumulative  Rod  Total  vs.  Time 


Figure  5.C.10.  “Pnp”  Configuration  1:  Instantaneous  Rod  Total  vs.  Time 
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nnn(Lab17)  area  Rod  Total  11788.4256  in  30000  sec 


Figure  5.C.1 1.  “nnn”  Configuration  1:  Cumulative  Rod  Total  vs.  Time 
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Figure  5.C.12.  “nnn”  Configuration  1:  Instantaneous  Rod  Total  vs.  Time 
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pnp(Lab17)  Rod  Total  vs.  Time 


Figure  5.C.14.  “pnp”  Configuration  1:  Instantaneous  Rod  Total  vs.  Time 
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Figure  5.C.15.  “nNp”  Configuration  1:  Cumulative  Rod  Total  vs.  Time 
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Figure  5.C.16.  “nNp”  Configuration  1:  Instantaneous  Rod  Total  vs.  Time 
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Figure  5.C.17.  “Ppp”  Configuration  1:  Cumulative  Rod  Total  vs.  Time 


Figure  5.C.18.  “Ppp”  Configuration  1:  Instantaneous  Rod  Total  vs.  Time 
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D.  CONFIGURATION  2  RESULTS 

The  combo  model’s  results  became  the  baseline  reference  for  all  data  comparison. 
The  power  required  to  achieve  steady  state  was  consistently  lower  than  any  other  model 
or  program,  regardless  of  initial  satellite  angular  rates.  The  advantage  of  lower  required 
power  levels,  however,  was  somewhat  diminished  by  the  apparent  trade-off  with 
increased  settling  times.  The  majority  of  the  other  models  and  programs  consistently 
performed  better  than  the  configuration  2  model  with  respect  to  time  to  achieve  steady 
state.  Nevertheless,  time  to  achieve  steady  state  has  been  considered  of  less  overall 
importance  than  power  use  for  the  initial  acquisition  of  orbit  attitude. 

All  32  initial  angular  rate  test  cases  were  analyzed  using  the  combo  model.  There 
was,  like  the  results  of  the  Bdot  program  runs,  no  apparent  correlation  between  sign  or 
magnitude  combinations  of  the  initial  angular  rates  and  the  resulting  time  or  power 
consumption  required  to  achieve  a  steady  state  condition.  There  appeared,  however,  to 
be  some  correlation  between  power  consumption  and  time  to  achieve  steady  state.  The 
initial  power  use  for  “high”  satellite  angular  rates  and  “low”  satellite  angular  rates  is 
comparable  between  test  cases,  regardless  of  which  angular  rate  had  the  highest 
magnitude.  All  configuration  2  test  cases  had  comparable  power  consumption  versus 
time,  once  the  initial  higher  rates  were  “knocked  down”  by  the  system.  Therefore,  the 
test  cases  that  took  the  longest  time  to  achieve  a  steady  state  condition,  correspondingly 
also  had  the  highest  power  consumption  levels.  Figures  5.D.1  through  5.D.18  show  the 
instantaneous  and  cumulative  power  consumption  results  for  configuration  2. 
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ppn(Combo)  area  Rod  Total  10532.4874  in  30000  sec 


Figure  5.D.I.  “ppn”  Configuration  2:  Cumulative  Rod  Total  vs.  Time 

ppn(Combo)  Rod  Total  vs.  Time 


Figure  5.D.2.  “ppn”  Configuration  2:  Instantaneous  Rod  Total  vs.  Time 
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Figure  5.D.3.  “NPn”  Configuration  2:  Cumulative  Rod  Total  vs.  Time 


Figure  5.D.4.  “NPn”  Configuration  2:  Instantaneous  Rod  Total  vs.  Time 
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NNp(Combo)  area  Rod  Total  57663.4151  ir  40000  sec 


Figure  5.D.5.  “NNp”  Configuration  2:  Cumulative  Rod  Total  vs.  Time 


NNp(Combo)  Rod  Total  vs.  Time 


Figure  5.D.6.  “NNp”  Configuration  2:  Instantaneous  Rod  Total  vs.  Time 
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pPn(Combo)  area  Rod  Total  77345.7856  in  80000  sec 


Time  (sec) 


Figure  5.D.7.  “pPn”  Configuration  2:  Cumulative  Rod  Total  vs.  Time 


Figure  5.D.8.  “pPn”  Configuration  2:  Instantaneous  Rod  Total  vs.  Time 


40 


Figure  5.D.9.  “Pnp”  Configuration  2:  Cumulative  Rod  Total  vs.  Time 

Pnp(Combo)  Rod  Total  vs.  Time 


Figure  5.D.10.  “Pnp”  Configuration  2:  Instantaneous  Rod  Total  vs.  Time 
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x  10 


nnn(Combo)  area  Rod  Total  65793.7608  in  120000  sec 


Figure  5.D.  1 1 .  “nnn”  Configuration  2:  Cumulative  Rod  Total  vs.  Time 
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Figure  5.D.12.  “nnn”  Configuration  2:  Instantaneous  Rod  Total  vs.  Time 
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pnp(Combo)  Rod  Total  vs.  Time 


Figure  5.D.14.  “pnp”  Configuration  2:  Instantaneous  Rod  Total  vs.  Time 
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Figure  5.D.15.  “nNp”  Configuration  2:  Cumulative  Rod  Total  vs.  Time 


Figure  5.D.16.  “nNp”  Configuration  2:  Instantaneous  Rod  Total  vs.  Time 
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Ppp(Combo)  area  Rod  Total  176876.4273  in  270000  sec 


Figure  5.D.17.  “Ppp”  Configuration  2:  Cumulative  Rod  Total  vs.  Time 

Ppp(Combo)  Rod  Total  vs.  Time 


Figure  5.D.18.  “Ppp”  Configuration  2:  Instantaneous  Rod  Total  vs.  Time 
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E.  CONFIGURATION  3  RESULTS 

When  the  Bdot  program  was  run  independently,  then  handed-off  to  the  LAB  17 
model,  using  the  final  conditions  from  the  Bdot  test  case  as  the  initial  conditions  for  the 
LAB  17  test  case,  the  following  results  occurred:  Compared  to  the  nine  combo  model  test 
cases  chosen  to  represent  the  full  spread  of  outcome  possibilities,  the  LAB  17  hand-off 
method  test  cases  tended  not  to  perform  as  well  compared  to  the  fastest  half  of  the  combo 
method  test  cases.  It  performed  better,  though,  than  the  slowest  half  of  the  combo 
method  test  cases.  Performance  levels,  mentioned  above,  include  both  power 
consumption  and  time  response  to  test  case  initial  conditions.  The  results  also  further 
amplified  the  finding  that  there  was  no  apparent  correlation  between  sign  or  magnitude  of 
the  test  case  initial  angular  rates  and  the  performance  of  the  system  in  power 
consumption  or  time  response.  Figures  5.E.1  through  5.E.2  show  the  instantaneous  and 
cumulative  power  consumption  results  for  configuration  3. 

Cases  “NPn,”  “NNp,”  “pPn,”  “nnn,”  and  “Ppp”  clearly  illustrate  the  power  spike 
occurring  at  the  hand-off  point  between  the  Bdot  and  Lab  17  models.  This  is  most  likely 
due  to  the  poor  initial  estimates  within  the  Kalman  Filter. 
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ppnhandoff(Lab17)  area  Rod  Total  86653.1201  in  90016  sec 


Figure  5.E.I.  “ppn”  Configuration  3:  Cumulative  Rod  Total  vs.  Time 
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Figure  5.E.2.  “ppn”  Configuration  3:  Instantaneous  Rod  Total  vs.  Time 
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Figure  5.E.3.  “NPn”  Configuration  3:  Cumulative  Rod  Total  vs.  Time 

NPnhandoff(Lab17)  Rod  Total  vs.  Time 


Figure  5.E.4.  “NPn”  Configuration  3:  Instantaneous  Rod  Total  vs.  Time 
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Figure  5.E.5.  “NNp”  Configuration  3:  Cumulative  Rod  Total  vs.  Time 


Figure  5.E.6.  “NNp”  Configuration  3:  Instantaneous  Rod  Total  vs.  Time 
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Figure  5.E.7.  “pPn”  Configuration  3:  Cumulative  Rod  Total  vs.  Time 


Figure  5.E.8.  “pPn”  Configuration  3:  Instantaneous  Rod  Total  vs.  Time 
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Figure  5.E.9.  “Pnp”  Configuration  3:  Cumulative  Rod  Total  vs.  Time 


Figure  5.E.  10.  “Pnp”  Configuration  3:  Instantaneous  Rod  Total  vs.  Time 
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nnnhandoff(Lab17)  area  Rod  Total  87620  848  in  110028  sec 


Figure  5.E.1 1.  “nnn”  Configuration  3:  Cumulative  Rod  Total  vs.  Time 


Figure  5.E.12.  “nnn”  Configuration  3:  Instantaneous  Rod  Total  vs.  Time 
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pnphandoff(Lab17)  area  Rod  Total  20813.7157  in  50010  sec 


Figure  5.E.13.  “pnp”  Configuration  3:  Cumulative  Rod  Total  vs.  Time 


Figure  5.E.14.  “pnp”  Configuration  3:  Instantaneous  Rod  Total  vs.  Time 
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nNphandoff(Lab17)  area  Rod  Total  213259.5242  in  69040  sec 


Figure  5.E.15.  “nNp”  Configuration  3:  Cumulative  Rod  Total  vs.  Time 


Figure  5.E.16.  “nNp”  Configuration  3:  Instantaneous  Rod  Total  vs.  Time 
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Figure  5.E.17.  “Ppp”  Configuration  3:  Cumulative  Rod  Total  vs.  Time 


Figure  5.E.18.  “Ppp”  Configuration  3:  Instantaneous  Rod  Total  vs.  Time 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

Table  6.A.1  summarizes  the  comparative  results  of  each  of  the  three  attitude 
control  methods  analyzed  in  this  study.  The  reference  combo  method  was  chosen 
primarily  for  its  relatively  superior  performance  in  power  consumption  required  to  reach 
a  steady  state.  The  nine  comparison  test  cases  were  chosen  to  represent  the  entire 
performance  range  of  the  combo  method  test  cases. 

The  most  obvious  conclusions  that  can  be  obtained  from  the  analysis  of  these 
control  methods  and  resulting  data  are  as  follows:  The  combo  model,  configuration  2 
generally  consumed  the  lowest  power  levels  of  any  other  method,  but  required  the 
longest  settling  times  for  medium  to  high  initial  angular  rates.  The  control  methods  that 
included  the  Bdot  control  law  (configuration  2  and  3)  quickly  reduced  the  satellite 
angular  rates  to  a  level  that  thence  required  only  minimal  power  consumption  to  achieve 
steady  state  pointing  accuracies  of  less  than  ten  degrees.  Time  to  achieve  system  steady 
state  was  achieved  more  quickly  and  with  less  power  consumption  for  cases  with  higher 
initial  angular  rates  when  the  hand-off  from  a  Bdot  control  law  occurred  after  a  threshold 
combined  angular  rate  was  achieved  instead  of  perfonning  the  hand-off  after  a 
predetermined  time  interval.  The  configuration  2  model  used  a  timed  hand-off  from  its 
internal  Bdot  law  control  loop,  and  the  configuration  3  model  handed-off  from  a  separate 
Bdot  program  to  the  various  other  models  after  a  threshold  combined  angular  rate  was 
achieved. 

The  first  of  the  three  main  analysis  objectives  was  to  examine  sensitivity  of  the 
three  configurations  to  initial  conditions.  All  three  configurations  perfonned  well  for  all 
the  test  case  initial  conditions  and  were  viable  options  for  an  ACS.  The  second  objective 
was  to  compare  the  differences  between  using  Bdot  or  no  Bdot  within  a  configuration. 
Although  there  was  questionable  improvement  in  system  time  responses  between  the 
configurations,  the  power  consumption  results  indicated  that  both  configurations  2  and  3, 
which  used  a  Bdot  control  loop,  required  less  power  to  achieve  steady  state  than 
configuration  1,  which  did  not  include  Bdot.  The  third  objective  was  to  explore  the 
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difference  between  using  a  timed  Bdot  hand-off  or  an  angular  rate  based  Bdot  hand-off. 
Configuration  2,  which  used  a  timed  hand-off,  had  better  power  consumption 
performance,  but  better  time  response  results  than  configuration  3,  which  used  the 
angular  rate  based  hand-off. 


Pwr  xl  e4 

Pwr  xl  e4 

Pwr  xl  e4 

Combo 

Time  req'd 

Hndff  bdot 

Lab17 

Init  Rates 

Labi  7 

ho  Lab17 

Combo 

Sim  sec 

Bdot  (sec) 

to  Lab17 

Simulink 

pPp 

20000 

9200 

m 

119278 

86653 

1 

30000 

20 

90000 

200000 

nnp 

30000 

3500 

NPp 

30000 

6700 

NPn 

169540 

162455 

49801 

30000 

6900 

56000 

67000 

nPn 

30000 

8900 

NNp 

5585031 

3814549 

57663 

40000 

11000 

190000 

92000 

Ppn 

50000 

3500 

Pnn 

60000 

1100 

Nnn 

70000 

3550 

Nnp 

70000 

3800 

pPn 

66757 

175825 

77346 

80000 

3500 

100000 

60000 

pNn 

80000 

8700 

NNn 

80000 

11000 

pnn 

90000 

380 

Npn 

90000 

840 

Pnp 

2474116 

103721 

77882 

100000 

650 

90000 

90000 

npn 

110000 

550 

npp 

110000 

640 

PPP 

110000 

3700 

NPP 

110000 

950 

PNn 

110000 

14000 

nnn 

1 

87621 

65794 

120000 

32 

110000 

30000 

PNp 

120000 

6700 

PNP 

120000 

8700 

PPp 

120000 

11000 

pnp 

98606 

20814 

84951 

140000 

12 

50000 

170000 

nPp 

150000 

8900 

nNp 

2677909 

2132595 

114607 

170000 

9044 

60000 

200000 

nNn 

170000 

9000 

PPn 

180000 

11200 

Ppp  1 

93255 

176876 

270000 

3000 

20000 

40000 

Table  6.A.  1 .  Test  Case  Table  Data. 


Figure  6.A.2  shows  common  system  power  consumption  response  characteristics 
that  can  be  gleaned  from  the  visual  comparison  of  the  three  configurations.  Magnetic 
torque  rod  power  consumption  was  lower,  and  high  angular  rates  were  generally  arrested 
faster,  when  a  Bdot  control  law  was  used  to  reduce  initial  angular  rates  to  a  more- 
manageable  “low”  level.  When  the  Bdot  control  law  was  used  to  achieve  a  sufficiently 
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low  threshold  combined  angular  rate,  as  opposed  to  a  timed  hand-off  from  a  Bdot  control 
law,  the  power  spike  that  often  occurs  at  the  hand-off  point  was  eliminated. 


Configuration  2 


Figure  6.A.2.  “pPn”  Test  Case. 


B.  RECOMMENDATIONS 

The  procedure  for  conducting  an  analysis  of  the  handoff  from  ‘nps.m,’  which  was 
the  Bdot  control  law  program  to  another  program  was  cumbersome  and  time-consuming. 
If  further  analysis  is  deemed  appropriate  to  determine  whether  a  Bdot  control  law  should 
be  used  for  NPSAT1,  a  combined  MATLAB  program  incorporating  each  of  the  smaller 
steps  should  be  developed  to  more  fully  automate  the  data  acquisition  and  analysis. 

It  appears,  howerer,  that  there  is  clear  evidence  that  the  Bdot  control  law  is  a 
powerful  and  efficient  method  to  arrest  initial  angular  rates,  regardless  of  their 
magnitude.  Should  it  be  decided  to  go  forward  with  some  implementation  of  a  Bdot 
control  law,  it  is  recommended  that  the  Bdot  control  law  be  a  selectable  control  loop 
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within  a  combined  program  or  system  model.  It  is  recommended  that  the  hand-off  that 
will  occur  from  the  Bdot  control  law  loop  should  be  done  after  a  threshold  combined 
angular  rate  is  achieved,  and  not  after  a  predetermined  time  interval.  The  time  for  the 
handoff  models  to  begin  steady  state  pointing  accuracy  acquisition  after  the  initial 
arresting  of  the  launch  vehicle  tip-off  rates  or  a  loss  of  attitude  control  can  be  minimized 
by  handing  off  control  as  soon  as  the  benefit  of  the  Bdot  control  law  ceases.  The  gains 
within  the  control  loop  that  will  receive  a  hand-off  from  the  Bdot  control  law  loop  can  be 
optimized  for  very  low  angular  rates,  since  the  Bdot  control  law  would  handle  any 
initially  higher  angular  rates. 

All  of  the  control  methods  analyzed  appeared  to  be  valid  control  methods  to 
achieve  three-axis  attitude  stabilization  using  only  magnetic  torquers  and  magnetometers 
for  active  control.  The  most  efficient  control  method  analyzed  appeared  to  be 
configuration  2.  This  method  was  achieved  by  simply  modifying  the  configuration  1 
model  to  perform  the  Bdot  hand-off  based  on  achievement  of  a  predetermined  combined 
angular  rate  as  opposed  to  a  predetermined  run  time.  This  method  also  appears  to 
eliminate  the  power  spike  that  seemingly  occurs  when  the  angular  rates  have  not  been 
sufficiently  arrested  prior  to  hand-off. 
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APPENDIX  A.  PROGRAMS 


NRL  BDOT  CONTROL  LAW  (NPS.M) 

%  Matlab  simulation  files  courtesy  of  Dr.  Glenn  Creamer,  NRL  Code  8231 
%  Main  program  (nps.m) 

clear  time  wl  w2  w3  ml  m2  m3  bl  b2  b3 

clear  uggl  ugg2  ugg3  umagl  umag2  umag3 

clear  yaw  angle  pitchangle  roll  angle  libration 

global  earth  radius  satellite  inertia  inv  satellite  inertia  dt 

%  output  storage  parameters 

count=0; 

k=0; 

store_count=3; 

%  conversion  factors 
d2r=pi/l  80; 
sqrt2=sqrt(2); 

%  Earth  parameters 
earth_radius=63  7 8.1363; 
earth_rate=7.272205217e-05; 
earth_gravity_constant=3 .9860044 1 5e+05 ; 

%  orbit  parameters 

inc=35*d2r;  %*****  UPDATE  DATA  ITEM 

ra=0*d2r; 

altitude=600;  %*****  UPDATE  DATA  ITEM 

radius=earth_radius+altitude ; 

%  orbit  rate  and  period 

orbit_rate=sqrt(earth_gravity_constant/radiusA3); 
peri  od=2  *  pi/orb  i  Irate; 

%  DCM  relating  perifocal  PQW  frame  to  Earth-centered  ECI  frame 
eci_from_pqw=[cos(ra)  -sin(ra)*cos(inc)  sin(ra)*sin(inc) 
sin(ra)  cos(ra)H5cos(inc)  -cos(ra)Hssin(inc) 

0  sin(inc)  cos(inc)]; 

%  initial  Greenwich  sidereal  time 
greenwich_time=0  *  d2r ; 

%  satellite  inertia  matrix  and  its  inverse 
satellite_inertia=[30  0  0  %*****  UPDATE  DATA  ITEM 
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031  0 

0  0  2]; 

inv_satellite_inertia=inv(satellite_inertia); 

%  satellite  torque  coil  capability 

m_max=30;  %*****  UPDATE  DATA  ITEM 

dipole_off=input('Disable  mag  torquers?  (l=Yes,0=No)'); 

%  magnetic  b-dot  control  feedback  gains 
k_mag=  1 00000000; 

%  initial  satellite  state 
true_anomaly=0  *d2r; 

%rates=[0.0 1  ;0.0 1  ;0.0 1  ]; 

rates=  [-0.1;  0.1;  0.01];  %  NPS  PDR  sim  results  used  these  rates 

%rates=  [  0.0;  0.001083;  0.0];  %  orbital  rate  @  600km  altitude  on  Y  axis 
%quats=[0;0;0;l]; 

quats=  [0.62721 137512625;  0.32650557562198;  %  ran  routine  initquat.m 

0.62721137512625;  0.32650557562198];  %  to  get  this 

state=[rates;quats] ; 
bfield_xyz_old=[0;0;0]; 

%  begin  time  loop 

dt=  2;  %*****  UPDATE  DATA  ITEM 

%t_final=3  600  *  .5 ; 

t_final=  input(’Input  sim  duration  (seconds):  ');  %  total  sim  time 
for  t=0:dt:t_final 

%  DCM  relating  spacecraft  XYZ  frame  to  ECI  frame 
xyz_from_eci=quats_to_dcm(quats); 

%  DCM  relating  rotating  orbit  frame 
%  (x  along  velocity,  z  along  zenith)  to  PQW  frame 

o rb_  1  ro m  p q w= [ - s i n ( truc  a n o m a  1  y )  cos(true  anomaly)  0 

0  0  1 
cos(true_anomaly)  sin(true_anomaly)  0]; 

%  DCM  relating  rotating  orbit  frame  to  ECI  frame 
orb_from_cci=orb_from_pqw*eci_from_pqw'; 

%  DCM  relating  spacecraft  XYZ  frame  to  rotating  orbit  frame 
xyz_from_orb=xyz_from_eciH5orb_from_eci'; 

%  3(yaw)-2(pitch)-l(roll)  Euler  angle  sequence  describing 
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%  spacecraft  frame  relative  to  orbit  frame 
yaw=atan(xyz_from_orb(  1 ,2)/xyz_from_orb(  1,1)); 
pitch=asin(-xyz_from_orb(  1 ,3)); 
roll=atan(xyz_from_orb(2 , 3  )/ xyz_from_orb(3 , 3  )) ; 

%  position  of  satellite  in  PQW  frame 
r _pqw=radius*  [eos(true_anomaly);sin(true_anomaly);0] ; 

%  position  of  satellite  in  ECI  frame 
r_eci=eci_from_pqw*r_pqw; 

%  position  of  satellite  in  ECEF  frame 
ecef_from_eci=[cos(greenwich_time)  sin(greenwich_time)  0 
-sin(greenwich_time)  cos(greenwich_time)  0 
0  0  1]; 
r_ecef=ecef_ffom_eci*r_eci; 

%  gravity  gradient  torque 

r_xyz=xyz_from_eci*r_eci; 

r_xyz_unit=r_xyz/sqrt(dot(r_xyz,r_xyz)); 

u_gg=3  *  orbit_rate  A2  *  skew(r_xyz_unit)  *  satellite_inertia  *  r_xyz_unit; 

%  angle  of  spacecraft  z-axis  relative  to  local  zenith  direction  (libration  angle) 
theta=acos(r_xyz_unit(3)); 

%  local  magnetic  field  in  ECEF  frame 
bfield_ecef=bfield(r_ecef,2); 

%  local  magnetic  field  in  XYZ  frame 
b  (i  c  1  d_x  y  z=x  y  z_  fro  m_e  c  i  *  cc  e  f_  fro  m_e  c  i '  *  b  (i  c  1  d_ec  c  f; 

%  time  derivative  of  local  magnetic  field  in  XYZ  frame 

bdot=(bfield_xyz-bfield_xyz_old)/dt; 

bfield_xyz_old=bfield_xyz; 

%  commanded  spacecraft  dipole 

if  dipole_off==  1 , 
dipole(l:3)=0; 
else 

ift  <  3600*30  %  may  want  to  remove  this  time-out 

dipole(l)=-k_mag*bdot(l); 
dipole(2)=-k_mag*bdot(2); 
dipole(3)=-k_mag*bdot(3); 
else 
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dipole(l:3)=0; 

end 

end 

%  dipole  saturation 
for  i=l:3 

if  abs(dipole(i))  >  mmax 
dipole(i)=m_max*sign(dipole(i)); 
end 
end 

%  corresponding  magnetic  torque 
u_mag=cross(dipole,bfield_xyz)'; 

%  store  critical  outputs  before  updating 
count=count+l; 
if  count  ==  store  count 
k=k+l; 
count=0; 
time(k)=t; 
wl(k)=rates(l); 
w2(k)=rates(2); 
w3(k)=rates(3); 
ml(k)=dipole(l); 
m2(k)=dipole(2); 
m3(k)=dipole(3); 
b  1  (k)=bfield_xyz(  1 ); 
b2(k)=bfield_xyz(2); 
b3(k)=bfield_xyz(3); 
uggl(k)=u_gg(l); 
ugg2(k)=u_gg(2); 
ugg3  (k)=u_gg(3 ) ; 
umag  1  (k)=u_mag(  1 ); 
umag2(k)=u_mag(2); 
umag3(k)=u_mag(3); 
yaw_angle(k)=yaw; 
pitch_angle(k)=pitch; 
roll_angle(k)=roll; 
libration(k)=theta; 
end 

%*****  PROPAGATE  ATTITUDE  AND  ORBIT  STATES  ***** 

%  updated  rotational  state  of  satellite  relative  to  ECI  frame 
u_ext=u_gg+u_mag ; 

state=satellite_dynamics(state,t,u_ext); 
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rates=state(  1:3); 
quats=state(4:7); 
quats=qnorm(quats); 

%  updated  greenwich  time 
greenwich_time=greenwich_time+earth_rate*dt; 

%  updated  true  anomaly 
true_anomaly=true_anomaly+orbit_rate*dt; 
end 

plotresult 

function  state=satellite_dynamics(state_old,t,u_ext) 
global  satellite  inertia  inv  satellite  inertia  dt 

rates_temp=state_old(  1:3); 

q_temp=state_old(4:7); 

rates_dot=inv_satellite_inertia* 

(u_ext-skew(rates_temp)Hssatellite_inertia*rates_temp); 
qmatrix=[q_temp(4) -q_temp(3)  q_temp(2)  q_temp(l); 

q_temp(3)  q_temp(4) -q_temp(l)  q_temp(2); 
-q_temp(2)  q_temp(l)  q_temp(4)  q_temp(3); 
-q_temp(l)  -q_temp(2)  -q_temp(3)  q_temp(4)]; 

quats_dot=0.5  *qmatrix*  [rates_temp;0] ; 
f  1  =dt*  [rates_dot;quats_dot] ; 

rates_temp=state_old(  1 :3)+f  1  ( 1 :3)/2; 
q_temp=state_old(4 : 7)+fl  (4 : 7)/2; 
rates_dot=inv_satellite_inertia* 

(u_ext-skew(rates_temp)H5satellite_inertia*rates_temp); 
qmatrix=[q_temp(4)  -q_temp(3)  q_temp(2)  q_temp(l); 

q_temp(3)  q_temp(4)  -q_temp(l)  q_temp(2); 
-q_temp(2)  q_temp(l)  q_temp(4)  q_temp(3); 
-q_temp(l)  -q_temp(2)  -q_temp(3)  q_temp(4)]; 

quats_dot=0.5  *qmatrix*  [rates_temp;0] ; 

G=dt*  [rates_dot;quats_dot] ; 

rates_temp=state_old(  1 :3)+f2(  1 :3)/2; 
q_temp=state_old(4 : 7)+12(4 : 7)/2; 
rates_dot=inv_satellite_inertiaH5 

(u_ext-skew(rates_temp)*satellite_inertiaH5rates_temp); 
qmatrix=[q_temp(4)  -q_temp(3)  q_temp(2)  q_temp(l); 

q_temp(3)  q_temp(4)  -q_temp(l)  q_temp(2); 
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-q_temp(2)  q_temp(l)  q_temp(4)  q_temp(3); 
-q_temp(l)  -q_temp(2)  -q_temp(3)  q_temp(4)]; 
quats_dot=0.5  *qmatrix*  [rates_temp;0] ; 
f3 =dt  *  [ratesdot;  quatsdot] ; 

rates_temp=state_old(  1 :3)+f3(l  :3); 
q_temp=state_old(4 : 7)+f3  (4 : 7) ; 

rates_dot=inv_satellite_inertia*(u_ext-skew(rates_temp)* 

satellite_inertia  *  rates_temp ); 
qmatrix=[q_temp(4)  -q_temp(3)  q_temp(2)  q_temp(l); 

q_temp(3)  q_temp(4)  -q_temp(l)  q_temp(2); 
-q_temp(2)  q_temp(l)  q_temp(4)  q_temp(3); 
-q_temp(l)  -q_temp(2)  -q_temp(3)  q_temp(4)]; 
quats_dot=0.5  *qmatrix*  [rates_temp;0] ; 
f4=dt*  [rates_dot;quats_dot] ; 

state=state_old+(f  1  +2  *  12+2  *  f3+f4)/ 6 ; 

function  bfield_ecef=bfield(r_ecef, order) 
global  earthradius 

if  order  >  8 
order=8; 

’Desired  order  exceeds  8  —  model  order  reset  to  8’ 
end 

%  current  orbit  position  in  Earth-fixed  REF  frame  (up,  south,  east) 
r=sqrt(r_ecef(l)A2+r_ecef(2)A2+r_ecef(3)A2); 
delta=asin(r_ecef(3)/r); 
cdelta=cos(delta); 
sdelta=sin(  delta); 
theta=pi/2-delta; 
stheta=sin(  theta); 
ctheta=cos(theta); 
phi=acos(r_ecef(  1 )/ (r  *  cdelta)); 
sphi=r_ecef(2)/ (r  *  c  delta) ; 
cphi=r_ecef(  1 )/ (r*  cdelta); 
if  sphi<0 
phi=-phi; 
end 

ref_from_ecef=  [cdelta*  cphi  cdelta*  sphi  sdelta 
sdelta*cphi  sdelta*  sphi  -cdelta 
-sphi  cphi  0]; 


%  S  factors 
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sO(l)=l; 

s(l,l)=l; 

if  order  >  1 
for  n=2:  order 
s0(n)=s0(n-l)*(2*n-l)/n; 
s(n,  1  )=sO(n)  *  sqrt(2  *n/(n+ 1 )); 
for  m=2:n 

s(nqn)=s(nqn-l)*sqrt((n-m+l)/(n+m)); 

end 

end 

end 

%  Legendre  functions  and  their  derivatives  wrt  theta 
pO(l)=ctheta; 
p(l,l)=stheta; 
dpO(l)=-stheta; 
dp(l,l)=ctheta; 
if  order  >  1 
for  n=2:  order 
if  n==2 

pO(2)=cthetaA2-l/3; 
dp0(2)=-2  *  stheta*  ctheta; 
else 

pO(n)=pO(n- 1 )  *  ctheta-(n- 1  )A2  /  ((2  *n- 1 )  *  (2  *n-3))  *p0(n-2); 
dpO(n)=dpO(n-l)*ctheta-pO(n-l)* 

stheta-(n- 1  )A2/((2  *n- 1  )*  (2  *n-3))  *  dp0(n-2) ; 
end 

p(n,n)=p(n- 1  ,n-  l)*stheta; 
dp(n,n)=dp(n- 1  ,n-  l)*stheta+p(n- 1  ,n- 1)*  ctheta; 
for  m=l:(n-l) 
if  m==(n-l) 

p(n,m)=p(n- 1  ,m)*ctheta; 
dp(n,m)=dp(n- 1  ,m)*ctheta-p(n- 1  ,m)*stheta; 
else 

p(n,m)=p(n- 1  ,m)*ctheta-((n-  l)A2-mA2)/((2*n-  l)*(2*n-3))*p(n-2,m); 
dp(n,m)=dp(n- 1  ,m)*ctheta-p(n- 1  ,m)* 

stheta-((n-l)A2-mA2)/((2*n-l)*(2*n-3))*dp(n-2,m); 

end 

end 

end 

end 

%  look-up  table  for  Gaussian  coefficients 

coeffs=gauss_coeffs; 

g=coeffs(  1 :44); 

h=coeffs(45:88); 
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%  first  order,  zeroth  degree  field  components 
br=2*(earth_radius/r)A3*sO(  l)*g(  l)*pO(  1); 
bt=-(earth_radius/r)A3  *  s0(  1 )  *  g(  1 )  *  dpO(  1 ); 
bp=0; 

%  higher  orders,  zeroth  degree  field  components 
if  order  >  1 

j=i; 

for  n=2:  order 

j=j+n; 

br=br+(n+l)*(earth_radius/r)A(n+2)*sO(n)*g(j)*pO(n); 

bt=bt-(earth_radius/r)A(n+2)*sO(n)*g(j)*dpO(n); 

end 

end 

%  all  orders,  higher  degree  field  components 

j=i; 

for  n=l:  order 

j=j+i; 

for  m=l:n 

br=br+(n+ 1 )  *  (earth_radius/ r)  A(n+2)  *  s(n,m)  *  (g(j )  * 

cos(m*phi)+h(j)*sin(m*phi))*p(n,m); 

bt=bt-(earth_radms/r)A(n+2)*s(n,m)*(g(j)*cos(m*phi)+h(j)* 

sin(m*phi))*dp(n,m); 
bp=bp-(  1  /  stheta)*  (earth_radius/r)A(n+2)  *m* 

s(n,m)*(-g(j)*sm(m*phi)+h(j)*cos(m*phi))*p(n,m); 

j=j+i; 

end 

end 

%  B-field  in  REF  frame 
bfield_ref=[br;bt;bp] ; 

%  B-field  in  ECEF  frame 
bfield_ecef=ref_from_ecef  *bfield_ref; 

function  coeffs=gauss_coeffs 

g=[-29682  -1789  -2197  3074  1685  1329  -2268  1249  769  941  782  291  -421 
1 16  -210  352  237  -122  -167  -26  66  64  65  -172  2  17  -94  78  -67 
1  29  4  8  10  -2  24  4  -1  -9  -14  4  5  0  -7]*1.0e-09; 
h=[0  5318  0  -2356  -425  0  -263  302  -406  0  262  -232  98  -301  0  44  157  -152 
-64  99  0  -16  77  67  -57  4  28  0  -77  -25  3  22  16  -23  -3  0  12  -20  7  -21 
12  10  -17  -10]*  1 .0e-09; 

coeffs=[g  h]; 
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function  qn=qnorm(q) 

mag=sqrt(dot(q,q)); 

qn=[q(l);q(2);q(3);q(4)]/mag; 


function  dcm=quats_to_dcm(q) 

dcm=[q(4)A2+q(l)A2-q(2)A2-q(3)A2  2*(q(l)*q(2)+q(3)*q(4)) 

2*(q(l)*q(3)-q(2)*q(4)); 

2*(q(l)*q(2)-q(3)*q(4))  q(4)A2-q(l)A2+q(2)A2-q(3)A2 

2*(q(2)*q(3)+q(l)*q(4)); 

2  *  (q(  1 )  *  q(3)+q(2)  *  q(4))  2*(q(2)*q(3)-q(l)*q(4)) 

q(4)A2-q(  1  )A2-q(2)A2+q(3)A2] ; 


function  skew_matrix=skew(vector) 
skew_matrix=[0  -vector(3)  vector(2); 

vector(3)  0  -vector(l); 
-vector(2)  vector(l)  0]; 


%plotresult.m 
r2d=  180/pi; 
figure 

subplot(31  l),plot(time, roll_angle*r2d), title(' Attitude  Error') 
ylabel('Roll(deg)'),grid 

subplot(312),plot(time,pitch_angle*r2d),ylabel('Pitch(deg)'),grid 
subplot(3 1 3),plot(time,yaw_angle*r2d),ylabel(’Y  aw(deg)’),grid 
xlabel(’time(sec)') 

figure 

subplot(31  l),plot(time,wl*r2d),title(’Body  Rates'),  ylabel(’wx(deg/sec)’),grid 
subplot(312),plot(time,w2*r2d),ylabel('wy(deg/sec)'),ylabel('wy(deg/sec)'),grid 
subplot(313),plot(time,w3*r2d),ylabel('wz(deg/sec)'),ylabel('wz(deg/sec)'),grid 
xlabel(’time(sec)') 


B.  NPS  (LAB15MOIDATA.M) 

%  Lab  15MOIdata.m  written  by  Prof  Barry  Leonard 

Ix=5;Iy=5.1;Iz=2;kax=.028;kay=.015;kaz=. 032*4;  %  save  for  b8 

%  steady  state  mag  act 

lxy=0  ;lxz=0  ;lyz=0 ; 
beta=0*pi/180; 

Imoi=[Ix  -Ixy  -Ixz;-Ixy  Iy  -Iyz;-Ixz  -Iyz  Iz]; 

Iinv=inv(Imoi); 
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pho=-.  1  ;tho=.  1  ;pso=-.  1 ; 
phdo=-0. 1  ;thdo=0. 1  ;psdo=0.0 1 ; 
sl=sin(pho/2);s2=sin(tho/2);s3=sin(pso/2); 
cl=cos(pho/2);c2=cos(tho/2);c3=cos(pso/2); 
elo=sl*c2*c3-cl*s2*s3;e2o=cl*s2*c3+sl*c2*s3;  %Wie pg.321 
e3o=cl*c2*s3-sl*s2*c3;e4o=cl*c2*c3+sl*s2*s3; 
Sl=sin(pho);S2=sin(tho);S3=sin(pso); 
Cl=cos(pho);C2=cos(tho);C3=cos(pso); 

wxo=phdo-psdo*S2-wo*S3*C2; 
wyo=thdo*Cl+psdo*C2*Sl-wo*(C3*Cl+S3*S2*Sl); 
wzo=psdo*C2*Cl-thdo*Sl-wo*(S3*S2*Cl-C3*Sl); 
eo=[elo  e2o  e3o  e4o];qo=eo; 

Ho=Imoi*[wxo  wyo  wzo]’; 

sat=30; 

gl=0.75;g2=0.82;g3=0.43;  %  varies  with  inclination 

%  (from  torque  ave.  with  high  sat.) 
ao=Ix-Iy+Iz;al=Iy-Iz;a2=Ix-Iz;a3=Iy-Ix; 


A=[0  0  0  1  0  0;0  0  0  0  1  0;0  0  0  0  0  1; 

-4*woA2*al/Ix  0  0  0  0  wo*ao/Ix; 

0  -3*woA2*a2/Iy  0  0  0  0; 

0  0  -woA2*a3/Iz  -wo*ao/Iz  0  0]; 

B=Kme*[0  0  0;0  0  0;0  0  0;kax/lx  0  0;0  kay/Iy  0;0  0  kaz/Iz]; 

Qx=l  10*eye(6);  Ru=l*[.l  0  0;0  .095  0;0  0  .008];%for  5  5.1  2  with  b8 
[K,S,e]=lqr(A,B,Qx,Ru); 

Aaa=A(  1:3,1 :3);Aab=A(  1 :3,4:6);Aba=A(4:6, 1 :3);Abb=A(4:6,4:6); 
Ba=B(  1:3,1 :3);Bb=B(4:6, 1 :3);Ka=K(  1:3,1 :3);Kb=K(  1 :3,4:6); 

Ga=G(  1:3,1 :3);Gb=G(4:6, 1 :3); 

Lr=l*diag([l  .5  .6]);  %for  5  5.1  2  withb8 

Lr=0.09*diag([l  .5  .6]); 

Ar=Abb-Lr*Aab-(Bb-Lr*Ba)*Kb; 

Br=Ar*Lr+Aba-Lr*Aaa-(Bb-Lr*Ba)*Ka; 

Cr=-Kb  ;Dr=-Ka-Kb  *  Lr; 

Kp=Ka;  Kd=Kb; 

sl=sin(pho/2);s2=sin(tho/2);s3=sin(pso/2); 
cl=cos(pho/2);c2=cos(tho/2);c3=cos(pso/2); 
elo=sl*c2*c3-cl*s2*s3;e2o=cl*s2*c3+sl*c2*s3;%Wie  pg.321 
e3o=cl*c2*s3-sl*s2*c3;e4o=cl*c2*c3+sl*s2*s3; 
Sl=sin(pho);S2=sin(tho);S3=sin(pso); 
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Cl=cos(pho);C2=cos(tho);C3=cos(pso); 

wxo=phdo-psdo  *  S2-wo  *83*02; 
wyo=thdo*Cl+psdo*C2*Sl-wo*(C3*Cl+S3*S2*Sl); 
wzo=psdo*C2*Cl-thdo*Sl-wo*(S3*S2*Cl-C3*Sl); 
eo=[elo  e2o  e3o  e4o];qo=eo; 

Ho=Imoi*[wxo  wyo  wzo]’; 

C.  NPS  (LAB16DATA.M) 

%Lab  1 6Data  written  by  Prof  Barry  Leonard 

clear 

Re=637 1 .2e3;mu=39860 1 .2e9;we=7.292  le-5;epsilon=l  1 .398*pi/l  80; 
Altitude=[450  500  550  600  650  700  750  800]*  le3;  %look  up  table  data  for  aero 
Density=[36.1  18  9.25  4.89  2.64  1.47  .837  ,439]*le-13;  %  lookup  table 

%  data  for  aero 

h=600e3;  %  spacecraft  altitude  (km) 

incln=(35)*pi/180;  %  inclination  (radians) 

beta=0*pi/180;  %  angle  between  outward  normal  vector  and  sun  vector 

nuo=-l  15*pi/l  80; 

alphao=0*pi/180; 

Lm=-70.454;Lgo=0;uo=(Lgo+Lm+90)*pi/180; 

wn=0;wo=sqrt(mu/(Re+h)A3); 

V=wo  *  (Re+h)  ;rho=asin(Re/(h+Re))  ;Kme=7 . 963  8e  1 5/(Re+h)A3 ; 
Cd=2.5;Kaero=0.5*Cd*VA2;  %  Cd  =  drag  coefficient 

psun=4 . 5  e-6  ;Psolar=2*psun; 

Area=[0.2674  0.2674  .1927];  dL=[.002  .002  .008]; 

lxy=0  ;lxz=0  ;lyz=0 ; 

Ix=5;Iy=5. 1  ;Iz=2;kax=l  ;kay=0. 1  *kax;kaz=kax*  1 ; 

Imoi=[Ix  -Ixy  -Ixz;  %  The  Inertia  tensor  or  Inertia  Matrix 
-Ixy  Iy  -Iyz; 

-Ixz  -Iyz  Iz] ; 

Iinv=inv(Imoi);  %  Inertia  Matrix  Inverse 

sat=30; 

gl=0.75;g2=0.82;g3=0.43;  %  varies  with  inclination 

%  (from  torque  ave.  with  high  sat.) 


ao=Ix-Iy+Iz;al=Iy-Iz;a2=Ix-Iz;a3=Iy-Ix; 

A=[0  0  0  1  0  0;0  0  0  0  1  0;0  0  0  0  0  l;-4*woA2*al/Ix  0  0  0  0  wo*ao/Ix;... 

0  -3*woA2*a2/Iy  0  0  0  0;0  0  -woA2*a3/Iz  -wo*ao/Iz  0  0]; 
B=Kme*[0  0  0;0  0  0;0  0  0;kax/lx  0  0;0  kay/Iy  0;0  0  kaz/Iz]; 

G=[0  0  0;0  0  0;0  0  0;1/Ix  0  0;0  1/Iy  0;0  0  1/Iz]; 

Qx=10*eye(6);  Ru=[.l  0  0;0  .1  0;0  0  .01]; 
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[K,S,e]=lqr(A,B,Qx,Ru); 

Aaa=A(  1:3,1 :3);Aab=A(  1 :3,4:6);Aba=A(4:6, 1 :3);Abb=A(4:6,4:6); 

Ba=B(  1:3,1 :3);Bb=B(4:6, 1 :3);Ka=K(  1:3,1 :3);Kb=K(  1 :3,4:6); 

Ga=G(  1:3,1 :3);Gb=G(4:6, 1 :3); 

Lr=0.2*eye(3);Ar=Abb-Lr*Aab-(Bb-Lr*Ba)*Kb; 

Br=Ar*Lr+Aba-Lr*Aaa-(Bb-Lr*Ba)*Ka; 

Cr=-Kb  ;Dr=-Ka-Kb  *  Lr; 

Kp=Ka;  Kd=Kb; 

pho=-.  1  ;tho=.  1  ;pso=-.0 1 ;  %  Euler  angles 

phdo=.  1  ;thdo=-0.0 1  ;psdo=0.0 1 ;  %  Euler  rates 

sl=sin(pho/2);s2=sin(tho/2);s3=sin(pso/2); 

cl=cos(pho/2);c2=cos(tho/2);c3=cos(pso/2); 

elo=sl*c2*c3-cl*s2*s3;e2o=cl*s2*c3+sl*c2*s3;  %  Wie  pg.321 

e3o=cl*c2*s3-sl*s2*c3;e4o=cl*c2*c3+sl*s2*s3; 

Sl=sin(pho);S2=sin(tho);S3=sin(pso); 

Cl=cos(pho);C2=cos(tho);C3=cos(pso); 

wxo=phdo-psdo*S2-wo*S3*C2;  %  wx 

wyo=thdo*Cl+psdo*C2*Sl-wo*(C3*Cl+S3*S2*Sl);  %  wy 

wzo=psdo*C2*Cl-thdo*Sl-wo*(S3*S2*Cl-C3*Sl);  %  wz 

eo=[elo  e2o  e3o  e4o]; 
qo=eo; 

Ho=Imoi*[wxo  wyo  wzo]’;  %  [Ho]  =  Angular  Momentum  [3x1] 

%SphericalHarmMag3 

r=Re+h;a=Re;  %  r  =  radius  from  sat  to  center  of  earth 

%  Re  =  a  =  Earth  radius  =  6371 .2e3  km 
kme=Kme*  le9;  %conversion  from  nT.mA3  to  T.mA3 

dth=5  ;tmin=0  ;tmax=  1 8  0 ;  Jmax=(tmax-tmin)/ dth; 
dphi=  1 0;phimin=- 1 80;phimax=l  80;Imax=(phimax-phimin)/dphi; 

ggg=[-30186  -2036  0  0  0  0  0  0  0;-1898  2997  1551  0  0  0  0  0  0;... 

1299  -2144  1296  805  0  0  0  0  0;951  807  462  -393  235  0  0  0  0;... 

-204  368  275  -20-161  -38  0  0  0;46  57  15  -210-1  -8  -114  0  0;... 

66  -57-7  7  -22  -9  11  -8  0;11  13  3  -12  -4  6  -2  9  1]; 

hhh=[0  5735  0  0  0  0  0  0  0;0  -2124  -37  0  0  0  0  0  0;... 

0  -361  249  -253  0  0  0  0  0;0  148  -264  37  -307  0  0  0  0;... 

0  39  142  -147  -99  74  0  0  0;0  -23  102  88  -43  -9  -4  0  0;... 

0  -68  -24  -4  1 1  27  -17  -14  0;0  4  -15  2  -19  1  18  -6  -19]; 
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gg=[-29615  -1728  0  0  0  0  0  0  0;-2267  3072  1672  0  0  0  0  0  0;... 

1341  -2290  1253  715  0  0  0  0  0;935  787  251  -405  110  0  0  0  0;. 
-217  351  222  -131  -169  -12  0  0  0;72  68  74  -161  -5  17  -91  0  0; 
79  -74  0  33  9  7  8  -2  0;25  6  -9  -8  -17  9  7  -8  -7]; 

hh=[0  5186  0  0  0  0  0  0  0;0  -2478  -458  0  0  0  0  0  0;... 

0  -227  296  -492  0  0  0  0  0;0  272  -232  1 19  -304  0  0  0  0;... 

0  44  172  -134  -40  107  0  0  0;0  -17  64  65  -61  1  44  0  0;... 

0  -65  -24  6  24  15  -25  -6  0;0  12  -22  8  -21  15  9  -16  -3]; 


for  J=l:Jmax+l; 

t=(J-l)*dth*pi/180;  th(J)=(J-l)*dth; 
Ct=cos(t);St=sin(t); 


for  I=l:Imax+l; 

p=(phimin+(I- 1 ) * dphi) *pi/ 180;  phi(I)=phimin+(I- 1 ) * dphi; 
Cp=cos(p);Sp=sin(p); 


for  i=l:8 
for  j=l:i+l 

Km(i,j)=((i- 1  )A2-(j  - 1  )A2)/ (2  *  i- 1  )/(2  *  i-3); 
if  i==  1 ;  Km=0;end 
end 
end 


Kin=tril(Kin,l); 
for  i=l:8 
for  j=l:i+l 
P(U)=Ct; 

if  i==l  &  j==2;P(i,j)=St*  1  ;end 
if  i>  1  &  j==i+ 1  ;P(i,j)=St*P(i- 1  ,j- 1  );end 
if  i==2;Pk=l;end; 
if  i>2;Pk=P(i-2,j);end 

if  i>l  &  j~=i+l;P(i,j)=Ct*P(i-l,j)-Km(i,j)*Pk;end 
end 
end 


P=tril(P,l); 
for  i=l:8 
for  j=l:i+l 
dp(U)=-St; 

if  i==l  &  j==2;dp(i,j)=St*0+Ct*  1  ;end 

if  i>  1  &  j==i+ 1  ;dp(i,j)=St*dp(i- 1  ,j- 1  )+Ct*P(i- 1  ,j- 1  );end 

if  i==2;dpk=0;end; 
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if  i>2;dpk=dp(i-2,j);end 

if  i>l  &  j~=i+l  ;dp(i,j)=Ct*dp(i-l  ,j)-St*P(i-l  j)-Km(i  j)*dpk;end 
end 
end 

dp=tril(dp,l); 
for  i=l:8 
for  j= 1:9 

dm=0;if  j==2;dm=  1  ;end 
if  j==l  &  i==l ;  S=l;end 
if  j=l  &  i~=  1  ;S(i,j)=S(i- 1  ,j)*(2*i- 1  )/i;end 
if  j>l;S(i,j)=S(i,j-l)*sqrt((i-(j-l)+l)*(dm+l)/(i+(j-l)));end 
end 
end 

S=tril(S,l); 

G=S.*gg; 

H=S.*hh; 
for  i=l:8 
j=l:i+l;m=j-l; 

F(i,j)=G(i,j).*cos(m*p)+H(i,j).*sin(m*p); 

end 


for  i=l:8 
j=l:i+l;m=j-l; 

f(i,j)=m.*(-G(i,j).*sin(m*p)+H(i,j).*cos(m*p)); 

end 

if  St==0;  St=0.001;  ;end 
for  i=l:8 

aa  1  (i)=(i+ 1  )*(a/r)A(i+2); 
aa2(i)=-aa  1  (i)/ (i+ 1 ); 
aa3(i)=aa2(i)/St; 
end 

FPS=sum((F.*P)')’  ;  Bl=aal’.*FPS;  Br8=sum(Bl);  Brl=Bl(l); 
FdpS=sum((F.!|sdp)')';  B2=aa2'.*FdpS;  Bt8=sum(B2);  Bt2=B2(l); 
fPS=sum((f.*P)')'  ;  B3=aa3’.*fPS;  Bp8=sum(B3);  Bp3=B3(l); 

br 8  (I , J)=Br 8/kme ;  br  1  (I , J)=Br  1  /kme ; 
bt8(I, J)=Bt8/kme;  bt  1  (I, J)=Bt2/kme; 
bp8(I, J)=Bp8/kme;  bp  1  (I, J)=Bp3/kme; 
b  1  (I,J)=(Br  1  ,A2+Bt2.A2+Bp3  A2)  A0. 5/kme; 
b8(I, J)=(Br8  A2+Bt8  A2+Bp8  A2)  A0. 5/kme; 
end 
end 
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Inclination=[0.01  10  20  30  40  50  60  70  80  90  100  1 10  120]*pi/180; 
gx=[.967  .955  .922  .876  .826  .78  .739  .709  .691  .686  .691  .709  .739]; 
gy=[0.0001  .0678  .256  .46  .632  .762  .857  .923  .965  .981  .965  .923  .857]; 
gz=[.804  .781  .711  .614  .522  .446  .39  .353  .335  .333  .335  .353  .39]; 


NPS  (LAB17.MDL) 

See  Figure  Al. 


Figure  Al.  Labl7.mdl 
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AA_BL1  mdl  To  Run:  open  &  run 
M-T>le  A_f«st_da1a1  m  men  A_second_data2.m 


E.  NPS  (NP SAT1ACSDAT A.M) 

%NPSAT1  ACSData  written  by  Prof  Barry  Leonard 

clear 

tic 

Re=6371.2e3;mu=398601.2e9;we=7.2921e-5;%earth  radius,  gravity  and  spin  rate 
Altitude=[450  500  550  600  650  700  750  800]*le3;%lookup  table  data  for  aero 
Density=[36.1  18  9.25  4.89  2.64  1.47  .837  .439]*le-13;%density  table  data 

h=600e3  ;incln=(35)*pi/ 1 80;beta=  1 5  *pi/ 1 80;%altitude,inclination,  solar  angle 
nuo=-l  15*pi/180;alphao=0*pi/180;%initial  subsolar  point  and  true  anomaly 
Lgo=0;%initial  position  of  Greenwich  Meridian  WRT  RAAN 
kpre=0;  %nodal  precession  constant  assumed  to  be  zero  here 
wn=kpre*(Re/(Re+h))A3.5*cos(incln);%nodal  precession  (zero  eccentricity) 
wo=sqrt(mu/(Re+h)A3);V=wo*(Re+h);  %orbital  angular  and  linear  velocity 
rho=asin(Re/(h+Re));  %earth  angular  radius 

Cd=2.5;psun=4.5e-6;%drag  coefficient  and  solar  pressure  constant~N/m/  2 
Kaero=0.5*Cd*VA2;Psolar=2*psun;  %  constants  for  aero  and  solar  torque  calc. 

Area=[0.2674  0.2674  .1927];  %projected  area~mA2  in  body  x,y,z  directions 
dL=[.002  .002  .008];  %predicted  (cp-cm)~m  in  body  x,y,z  directions 

Ix=5;Iy=5.1;Iz=2;Ixy=0;Ixz=0;Iyz=0;  %  moments  of  inertia  (MOI)~kg.mA2 
Imoi=[Ix  -Ixy  -Ixz;-Ixy  Iy  -Iyz;-Ixz  -Iyz  Iz];  %MOI  matrix 
Iinv=inv(Imoi);  %MOI  matrix  inverse 

ao=Ix-Iy+Iz;al=Iy-Iz;a2=Ix-Iz;a3=Iy-Ix;%MOI  combination  constants 

pho=-.l;tho=.l;pso=-0.1;  %initial  Euler  angles  (phi, theta, psi)  ~r 
phdo=-0. 1  ;thdo=0. 1  ;psdo=0.0 1 ;  %  initial  Euler  angle  rates  ~r/s 

%calculation  of  initial  quaternion  (qo)  and  angular  momentum  (Ho): 
sl=sin(pho/2);s2=sin(tho/2);s3=sin(pso/2);cl=cos(pho/2);c2=cos(tho/2); 
c3=cos(pso/2);  qlo=sl  *c2*c3-c  1  *s2*s3;q2o=c  1  *s2*c3+sl  *c2*s3;%Wie  pg.32 1 
q3o=cl*c2*s3-sl*s2*c3;q4o=cl*c2*c3+sl*s2*s3; 

Sl=sin(pho);S2=sin(tho);S3=sin(pso);Cl=cos(pho);C2=cos(tho);C3=cos(pso); 
wxo=phdo-psdo*S2-wo*S3*C2; 
wyo=thdo*Cl+psdo*C2*Sl-wo*(C3*Cl+S3*S2*Sl); 
wzo=psdo*C2*Cl-thdo*Sl-wo*(S3*S2*Cl-C3*Sl); 
qo=[qlo  q2o  q3o  q4o]; 

Ho=Imoi*[wxo  wyo  wzo]’; 

%Calculation  of  spherical  harmonic  magnetic  field  model  (Wertz  pp779-783) 
r=Re+h;a=Re;  %definitions  to  match  Wertz  pg  780 
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%  IGRF  Epoch  2000  Guassian  coefficients  ~  n.T  : 
gg=[-29615  -1728  0  0  0  0  0  0  0;-2267  3072  1672  0  0  0  0  0  0;... 

1341  -2290  1253  715  00  00  0;935  787  251  -405  110  00  00;... 

-217  351  222  -131  -169  -12  0  0  0;72  68  74  -161  -5  17  -91  0  0;... 

79  -74  0  33  9  7  8  -2  0;25  6  -9  -8  -17  9  7  -8  -7]; 
hh=[0  5186  0  0  0  0  0  0  0;0  -2478  -458  0  0  0  0  0  0;... 

0  -227  296  -492  0  0  0  0  0;0  272  -232  1 19  -304  0  0  0  0;... 

0  44  172  -134  -40  107  0  0  0;0  -17  64  65  -61  1  44  0  0;... 

0  -65  -24  6  24  15  -25  -6  0;0  12  -22  8  -21  15  9  -16  -3]; 

Kme=(a/r)A3  *sqrt(gg(  1 ,  l)A2+gg(  1 ,2)A2+hh(  1 ,2)A2)*  1  e-9;%dipole  strength-Tesla 

ka=[. 02863  .01534  . 1 309] ;sat=30;%actuator  gains  and  saturation 
%ka=[.03  .015  .13]; 

%Calc  of  state  space  A,B,G  matrices  follows  (xdot=Ax+Bu+Gw,w=disturbance) 
A=[0  0  0  1  0  0;0  0  0  0  1  0;0  0  0  0  0  l;-4*woA2*al/Ix  0  0  0  0  wo*ao/Ix;... 

0  -3*woA2*a2/Iy  0  0  0  0;0  0  -woA2*a3/Iz  -wo*ao/Iz  0  0]; 

B=Kme*[0  0  0;0  0  0;0  0  0;ka(l)/lx  0  0;0  ka(2)/Iy  0;0  0  ka(3)/Iz]; 

G=[0  0  0;0  0  0;0  0  0;1/Ix  0  0;0  1/Iy  0;0  0  1/Iz]; 

Qx=l  10*eye(6);  Ru=l*[.l  0  0;0  .095  0;0  0  .008];  %LQR  weighting  matrices 
[K,S,e]=lqr(A,B,Qx,Ru);  %LQR  gain  calculation 

%partitioning  of  A,B,G  matrices  required  for  reduced  order  estimator; 

Aaa=A(  1:3,1 :3);Aab=A(  1 :3,4;6);Aba=A(4;6, 1 :3);Abb=A(4:6,4:6); 

Ba=B(l;3,l  :3);Bb=B(4;6, 1 :3);Ka=K(  1 :3, 1 :3);Kb=K(l  ;3,4;6); 

Ga=G(  1:3,1 :3);Gb=G(4:6, 1 :3); 

Lr=0.09*diag([l  .5  .6]);  %Estimator  gain  found  by  simulation 

del_t=5;tmin=0;tmax=180;%colatitude  increment,  min,  max  (theta) 
Jmax=(tmax-tmin)/del_t;  %  max  value  of  J  index 
del_p=10;pmin=-180;pmax=180;  %longitude  increment,  min,  max  (phi) 
Imax=(pmax-pmin)/del_p;  %max  value  of  I  index 

for  J=l:Jmax+l; 

t=(J-l)*del_t*pi/180;Ct=cos(t);St=sin(t);  %sine  &  cosine  of  theta 
th(J)=(J-l)*del_t;  %theta  vector  (look  up  tables  Mag  Fid  Model) 
for  I=l:Imax+l; 

p=(pmin+(I-l)*del_p)*pi/180;Cp=cos(p);Sp=sin(p);%sin  &  cos  phi 
phi(I)=pmin+(I-l)*del_p;%phi  vector  (look  up  tables  Mag  Fid  Model) 

%  Calculation  of  Legendre  function  constant  (Km)  follows 
for  i=l:8 
for  j=l:i+l 

Km(i,j)=((i- 1  )A2-(j  - 1  )A2)/ (2  *  i- 1  )/(2  *  i-3) ; 
if  i==  1 ;  Km=0;end 
end 
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end 

Kin=tril(Km,l); 


%Calculation  of  Legendre  function  (P)  follows 
for  i=l:8 
for  j=l:i+l 
P(U)=Ct; 

if  i==l  &  j==2;P(i,j)=St*  1  ;end 

if  i>  1  &  j==i+ 1  ;P(i,j)=St*P(i- 1  ,j- 1  );end 

if  i==2;Pk=l;end; 
if  i>2;Pk=P(i-2,j);end 

if  i>l  &  j~=i+l;P(i,j)=Ct*P(i-l,j)-Km(i,j)*Pk;end 
end 
end 

P=tril(P,l); 

%  Calculation  of  partial  of  Legendre  function  WRT  theta  (dp)  follows 
for  i=l:8 
for  j=l:i+l 
dp(U)=-St; 

if  i==l  &  j==2;dp(i,j)=St*0+Ct*  1  ;end 

if  i>  1  &  j==i+ 1  ;dp(i,j)=St*  dp(i- 1  ,j- 1  )+Ct*P(i- 1  ,j- 1  );end 

if  i==2;dpk=0;end; 

if  i>2;dpk=dp(i-2,j);end 

if  i>l  &  j~=i+l  ;dp(i,j)=Ct*dp(i-l  ,j)-St*P(i-l  ,j)-Km(i,j)*dpk;end 
end 
end 

dp=tril(dp,l); 

%  Calculation  of  Gausian  coefficient  nonn.  factor  (S)  follows 
for  i=l:8 
for  j=  1:9 

dm=0;if  j==2;dm=  1  ;end 
if  j==  1  &  i==l ;  S=l;end 
if  j==l  &  i~=  1  ;S(i,j)=S(i- 1  ,j)*(2*i- 1  )/i;end 
if  j>l;S(i,j)=S(i,j-l)*sqrt((i-(j-l)+l)*(dm+l)/(i+(j-l»);end 
end 
end 

S=tril(S,l); 

%Calculation  of  normalized  Gaussian  coefficients  ~  n.T  : 

Gg=S.*gg;  Hh=S.*hh; 


for  i=l:8 
j=l:i+l;m=j-l; 
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F(i,j)=Gg(i,j).H5cos(m*p)+Hh(i,j).*sin(m*p);  %definition 
end 

for  i=l:8 
j=l:i+l;m=j-l; 

f(i,j)=m.*(-Gg(i,j).*sin(m*p)+Hh(i,j).*cos(m*p));%  definition 
end 

if  St=0;  St=0.001;  ;end  %singularity  avoidance  in  following  loop 
%Calculation  of  altitude  dependent  coefficients  follows 
for  i=l:8 

aal(i)=(i+l)*(a/r)A(i+2);  aa2(i)=-aal(i)/(i+l);  aa3(i)=aa2(i)/St; 
end 

%Calculation  of  radial, theta  and  phi  magnetic  field  components  -Tesla 
FPS=sum((F.*P)')'  ;  Bl=aal'.*FPS;  Br8=sum(Bl)*le-9; 
FdpS=sum((F.*dp)')';  B2=aa2’.*FdpS;  Bt8=sum(B2)*le-9; 
fPS=sum((f.*P)')'  ;  B3=aa3’.*fPS;  Bp8=sum(B3)*le-9; 

%Calculation  of  normalized  field  components  (input  to  lookup  tables) 

br8(I,J)=Br8/Kme; 

bt8(I,J)=Bt8/Kme; 

bp8(I,J)=Bp8/Kme; 

%Calculation  of  field  vector  normalized  magnitude: 
b8(I,  J)=(Br8  ,A2+Bt8 .  A2+Bp8  ,A2).A0.5/Kme; 

end  %end  I  loop 
end  %end  J  loop 

%  Table  inputs  for  average  actuator  gains  vs  inclination  follows: 
Inclination=[0.01  10  20  30  40  50  60  70  80  90  100  1 10  120]*pi/180; 
gx=[.967  .955  .922  .876  .826  .78  .739  .709  .691  .686  .691  .709  .739]; 
gy=[0.0001  .0678  .256  .46  .632  .762  .857  .923  .965  .981  .965  .923  .857]; 
gz=[.804  .781  .711  .614  .522  .446  .39  .353  .335  .333  .335  .353  .39]; 
toe 

F.  NPS  (NPSATACSINTRVW.MDL) 

See  Figure  A2. 
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Figure  A2.  NPSATACSIntRvw.mdl. 
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G.  NPS  (GRABDATASET.M) 

%grabdataset.m  written  by  Lt  Todd  Zirkle  for  manipulating  data  captured  from 
%  a  SIMULINK  run  during  analysis  for  NPSSAT 1  NPS  Monterey 

% 

%  This  program  measures  and  plots  the  instantaneous  and  cumulative  total  of 
%  commanded  rod  torques  and  the  time  required  to  achieve  a  3 -axis  steady 
%  state  pointing  accuracy  of  less  than  10  degrees. 

%  This  program  also  solves  and  plots  the  area  under  the  instantaneous  rod 
%  total  curve  to  produce  an  energy  index  for  comparison  with  other  runs 
%  using  different  initial  conditions. 

% 

%  Eleven  variables  are  recorded  in  a  matrix  as  a  ‘.mat  file’  from 
%  a  simulation  using  one  of  two  modified  SIMULINK  simulations  originally 
%  written  by  Prof  Barry  Leonard. 

%  ‘Labl7.mdl’  05/18/01  and  ‘NPSATACSIntRvw.mdT  07/16/01 

% 

%  This  program  picks  off  some  of  the  desired  variables. 

%  The  desired  torque  rod  outputs  in  AmA2  are  recorded  as 
%  ‘rodl’,  ‘rod2’,  and  ‘rod3’ 

%  The  total  commanded  torque  from  the  three  rods  combined  is  ‘rodtot’ 

% 

%  After  loading  the  desired  ‘.mat’  file  containing  the  data  matrix 
%  into  the  MATLAB  workspace,  the  matrix  size  is  used  to  edit  this 
%  program’s  array  sizes. 

clear  ans  xacksis  DataArray  rodone  rodl  rod2  rod3  rodtot 

clear  rodtwo  rodthree  rodtotal 

clear  PhdT  ThdT  PsdT  PhT  ThT  PsT 

clear  Phd  Thd  Psd  Ph  Th  Ps  lasttime 

clear  Timegap  RulerTime  miniarea  areatot  areaCumPlot 

DataArray  =  NPn7600';  %  matrix  name  from  workspace  containing  data 

xacksis  =  DataArray(l:11370,l);  %  array  size  must  be  updated  to 

%  match  matrix  size 

xaxiss=xacksis'; 

rodl=  DataArray(l:  11370,2);  rodone  =  rodl'; 
rod2=  DataArray(  1 : 1 1370,3);  rodtwo  =  rod2’; 
rod3=  DataArray(  1 : 1 1370,4);  rodthree  =  rod3’; 
rodtot=  DataArray(l :  1 1370,5);  rodtotal  =  rodtot’; 

PhdT=  DataArray(l:  1 1370,6);  Phd  =  PhdT’; 

ThdT=  DataArray(  1 : 1 1370,7);  Thd  =  ThdT’; 

PsdT=  DataArray(l:  11370,8);  Psd  =  PsdT’; 

PhT=  DataArray(l:  1 1370,9);  Ph  =  PhT’; 

ThT=  DataArray(l :  1 1370,10);  Th  =  ThT’; 

PsT=  DataArray(l:  11370,11);  Ps  =  PsT’; 
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areatot=0;  %  essentially  starts  the  areatot  at  the  origin  which  is  zero 

for  k  =  1 : 1 1 369  %  will  be  one  less  than  the  size  of  the  matrix  dataset 
Timegap(k)=  xacksis(k+l)-xacksis(k);  %  delta  t  for  integrating  area  under  curve 
holdtotal=areatot; 

miniarea(k)=rodtot(k)*Timegap(k);  %  represents  area  under  curve 

%  from  (t+n)  to  (t+n+timegap) 
areatot=holdtotal+miniarea(k);  %  cumulative  area  up  til  this  point 
areaCumPlot(k)=areatot;  %  records  cumulative  area  at  each  time  interval 

end 

areaCumPlot(  1 1370)=areaCumPlot(l  1369);  %  the  last  delta  t  is  reused  for 

%  the  last  miniarea  calculation 

lasttime=xacksis(  1 1370);  %  last  time  of  simulation  recorded  for  plot  label  use 

figure  %  plots  rod  total  commanded  torques  vs.  time 

plot(xaxiss,rodtot,'r'),title(’NPn  (Labl7)  Rod  Total  vs.  Time'), 
ylabel(’(AmA2)'),grid 
xlabel(’Time  (sec)') 

figure  %  plots  cum  area  under  curve  vs.  time  [energy  index  (AmA2  *  sec)] 

plot(xaxiss,areaCumPlot,'r'), 

title(['NPn(Labl7)  area  Rod  Total  ’,num2str(areatot),’  in  ’,num2str(lasttime),'  sec’]) 
ylabel('  '),grid 
xlabel(’Time  (sec)') 


H.  NPS  (CREATE  VARIABLES.M) 

%create  variables  written  by  Lt  Todd  Zirkle  for  capturing  data  from 
%  NRL  program  written  by  Dr.  Glenn  Creamer,  NRL  Code  823 1 

%  (nps.m) 

%  This  program  captures  1 1  variables  for  analysis  and  plotting  purposes 
%  Before  running  this  program,  run  ‘nps.m’  after  ensuring  the  initial  conditions 
%  have  been  set  to  the  desired  values. 

%  This  program  will  capture  the  variables  by  assigning  them 
%  to  new  variable  names. 

%  After  renaming  the  1 1  variables,  the  program  clears  all  ‘nps.m’  variables. 

%  The  resulting  data  is  saved  by  saving  the  MATLAB  workspace 
%  to  a  data  folder. 

%  It  is  recommended  that  the  ‘.mat’  file’s  name  in  some  way  reflects 
%  the  initial  conditions  used  to  produce  the  data,  since  the  only  other 
%  way  to  ascertain  the  initial  conditions  is  to  pick  out  the  first  row  of 
%  data  from  the  six  variables  assigned  to  rates  and  angles. 

%  EXAMPLE  FILE  NAME:  ‘bdot_pNp_testcasel.mat’ 

%  In  this  example,  pNp  stands  for 

%  ‘phidot’  =.01  rad/s,  ‘thdot’  =-.l  rad/s,  ‘psidot’  =.01  rad/s 
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clear  brodl  brod2  brod3  brodtot  bwwl  bww2  bww3  bpitt  brole  byeaw 

brodl=[ml];  %  command  to  torque  rod  1  (phi)  ((from  0  to  30  AmA2)) 

brod2=[m2];  %  command  to  torque  rod  2  (theta)  ((from  0  to  30  AmA2)) 

brod3=[m3];  %  command  to  torque  rod  3  (psi)  ((from  0  to  30  AmA2)) 

brodtot=[rod3tot];  %  rod3tot  =  the  sum  of  [abs(rodl)+abs(rod2)+abs(rod3)] 

bwwl=[wl]; 

bww2=[w2]; 

bww3=[w3]; 

bpitt=[pitch_angle]; 

brole=  [rollangle] ; 

byeaw=[yaw_angle]; 

btime=[time]; 

clear  yaw  xyzfromorb  xyzfromeci  u  ext  u_gg  umag 

clear  watt  hrs  total  t  t_  final  theta  trueanomaly 

clear  store  count  state  sqrt2  satellite  inertia  inv  satellite  inertia 

clear  ra  radius  rates  roll  r  eci  r_pqw  r  xyz  r  xyz  unit 

clear  pitch  quats  r2d  r  ecef  orbfromeci  orbfrompqw 

clear  orbit  rate  period  inv  satellite  k  kmag  mmax 

clear  eci_from_pqw  greenwich_time  i  inc  earth_gravity_constant 

clear  earthradius  earth  rate  eceffrom  eci  d2r  dipole 

clear  dipole  off  dt  count  bdot  bfield  ecef  bfield  xyz 

clear  bllchl  xyz  old  ans  ampsec  total  altitude 

clear  ecef  from  eci 

clear  time  wl  w2  w3  w4  ml  m2  m3  bl  b2  b3 
clear  rate_goal  rg  time_counter  pwrl  pwr2  pwr3  pwr4  pwrSS 
clear  pwr6  pwr7  pwr8  percent_pwr  rod_current  rod3tot 
clear  pwr_u  uggl  ugg2  ugg3  umagl  umag2  umag3 
clear  yaw  angle  pitch  angle  roll  angle  libration 


I.  NPS  (NPSSAT1TESTCASE.M) 

%NPSsatl  Test  case  data  manipulator 

% 

%  Written  by  Lt  Todd  Zirkle  for  manipulating  data  captured  from 
%  MATLAB  programs  and  SIMULINK  runs  during  analysis 
%  for  NPS S ATI  NPS  Monterey 

% 

%  This  program  computes  and  plots  the  results  of  combining 
%  two  different  simulations’  data.  A  Bdot  control  law  simulation  is 
%  used  to  arrest  rotation  rates  up  to  a  predetermined  rate  goal,  then  the 
%  final  states  of  that  simulation  are  used  for  the  initial  conditions  of  the 
%  second  program,  which  further  arrests  the  rotation  rates  and  then 
%  achieves  a  predetennined  three-axis  pointing  accuracy. 
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% 

%  Uses  an  NRL  program  written  by  Dr.  Glenn  Creamer,  NRL  (nps.m) 

%  Uses  two  modified  SIMULINK  simulations  originally  written 
%  by  Prof  Barry  Leonard.  ‘Labl7.mdl’  05/18/01 

%  ‘NPSATACSIntRvw.mdl’  07/16/01 

% 

%  This  program  essentially  combines  two  other  programs  written  by  Lt  Zirkle: 
%  ‘grabdataset.m’  and  ‘create_variables.m’ 

%  As  the  original  programs  gathered  rodtotal  and  energy  index  vs.  time  plots 
%  for  the  Dr.  Creamer  program  and  the  Prof  Leonard  programs,  this  program 
%  does  the  same  for  the  case  of  both  programs  being  run  in  a  handoff  fashion 
%  as  described  above. 

% 

%  This  program  captures  torque  rod  total  and  time  variables  for  analysis  and 
%  plotting  purposes 

%  Before  running  this  program,  run  ‘nps.m’  after  ensuring  the  initial  conditions 
%  have  been  set  to  the  desired  values. 

%  This  program  will  capture  the  variables  by  assigning  them  to  new  variable 
%  names.  After  renaming  the  variables,  the  program  clears  all  non-essential 
%  variables. 

% 

%  This  program  measures  and  plots  the  instantaneous  and  cumulative  total  of 
%  commanded  rod  torques  and  the  time  required  to  achieve  a  3 -axis  steady 
%  state  pointing  accuracy  of  less  than  10  degrees. 

%  This  program  also  solves  and  plots  the  area  under  the  instantaneous  rod  total 
%  curve  to  produce  an  energy  index  for  comparison  with  other  runs  using 
%  different  initial  conditions. 

% 

%  After  loading  the  desired  ‘.mat’  files  containing  data  matrices  into  the 
%  MATLAB  workspace,  the  matrix  sizes  are  used  to  edit  this  program’s  array 
%  sizes. 

%  The  data  are  saved  as  a  ‘.mat’  file.  One  of  the  two  Prof  Leonard  programs 
%  is  run  after  initially  running  the  Bdot  program.  This  data  is  stored  as 
%  “handoff  data” 

%  first  import  handoff  data,  then  bdot  data 
%  ensure  handoff  variable  name  appears  on  the 
%  line  that  creates  the  HDataArray  matrix 
%  then  ensure  the  matrix  sizes  are  updated 

clear  Hxacksis  HDataArray  Hrodone  Hrodl  Hrod2  Hrod3  Hrodtot 

clear  Hrodtwo  Hrodthree  Hrodtotal 

clear  HPhdT  HThdT  HPsdT  HPhT  HThT  HPsT 

clear  HPhd  HThd  HPsd  HPh  HTh  HPs 

clear  HTimegap  HRulerTime  Hminiarea  Hareatot  HareaCumPlot 
clear  newtime  CompositeRodtotal  Composite  Time  bdottime 
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HDataArray  =  NPnhandoff6600’;  %  input  name  of  handoff  matrix 
Hxaxis  =  HDataArray(  1:8135,1);  %  input  size  of  handoff  matrix 

%  Hxaxis  will  be  the  (time  (x  axis)) 

Hxacksis  =  Hxaxis'; 

bdottime  =  btime';  %  time  (x-axis)  before  handoff 

Hrodtot=  HDataArray(l:8135,5);  Hrodtotal  =  Hrodtot'; 

HTimegap(  1  )=  Hxacksis(2)-Hxacksis(  1 ); 

newtime(l)=6908;  %  final  bdot  time  must  be  entered  here 

for  k=l :  1 1 5 1  %  input  matrix  size  of  bdot  matrix 

Composite_Rodtotal(k)=brodtot(k); 

Composite_Time(k)=bdottime(k); 

end 

Composite_Rodtotal(  1 1 52)=Hrodtot(  1 ); 

Composite_Time(l  152)=newtime(l); 

for  k  =  2:8 134  %  will  be  one  less  than  handoff  matrix  size 

HT  imegap(k)=  Hxacksis(k+ 1  )-Hxacksis(k); 
newtime(k)=newtime(k- 1  )+HT  imegap(k); 

Composite_Rodtotal(  115  l+k)=Hrodtot(k);  %  input  one  less  than  bdot  size 
Composite_Time(  115  l+k)=newtime(k);  %  input  one  less  than  bdot  size 
end 

HTimegap(8135)=  Hxacksis(8135)-Hxacksis(8134);  %  h-off  x2,  h-off  less  one 
newtime(8135)=newtime(8134)+HTimegap(8135);  %  h-off,  h-off  less  one,  h-off 
Composite_Rodtotal(l  151+8135)=Hrodtot(8135);  %  bdot  plus  h-off  size,  h-off 

Composite_Time(l  151+8135)=newtime(8135);  %  bdot  plus  h-off  size,  h-off 

figure  %  plots  composite  rod  total  vs.  time 

plot(Comp_Time,Comp_Rodtotal,'r'),title('Handoff  Rod  Total  vs  Time'), 
ylabel(’  AmA2  '),grid 
xlabel(’Time  (sec)') 


J.  NPS  (RODTOTAL  VS_TIME.M) 

%Rodtotal_vs_Time  written  by  Lt  Todd  Zirkle  to  plot  simulated  torque  rod  data 
%  for  analysis  of  simulations  supporting  NPSSAT1  NPS  Monterey 

% 

%  This  program  captures  a  row  of  data  from  a  saved  1 1  row  matrix 
%  The  matrix  was  produced  from  a  simulation  using  one  of  two  modified 
%  SIMULINK  simulations  originally  written  by  Prof  Barry  Leonard. 

%  ‘Labl7.mdT  05/18/01  and  ‘NPSATACSIntRvw.mdT  07/16/01 
%  This  program  could  be  used  to  capture  any  row  for  plotting  purposes 
%  This  program  grabs  torque  rod  totals  and  plots  them  vs.  time 
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clear  ans  HDataArray  Hxaxis  Hxacksis  Hrodtot 


HDataArray  =  NPnCombo6908’;  %  input  matrix  name  from  ‘.mat’  file 

%  loaded  into  workspace  prior  to  run 

Hxaxis  =  HDataArray(  1:1720,1);  %  input  array  size  here  from  matrix  size 
Hxacksis  =  Hxaxis';  %  time  axis 

Hrodtot=  HDataArray(  1:1720,5);  Hrodtotal  =  Hrodtot'; 

figure  %  plots  torque  rod  total  (sum  of  absolute  values  of  3  rods)  vs  time 

plot(Hxaxis, Hrodtot, 'r'),title(’B. Leonard  Combo  with  Bdot  Rod  Total  vs  Time'), 
ylabelC  AmA2  '),grid 
xlabel(’Time  (sec)') 

K.  NPS  (PLOTRESULT2.M) 

%plotresult2.m  written  by  Todd  A.  Zirkle  for  graphical  analysis 
%  of  NPSSAT1  data  NPS  Monterey 


figure 

subplot(41  l),plot(time,wl), 

title(['w_t_o_t_a_l  <  ’,num2str(rate_goal),'  degrees  (’,num2str(rg),'  rad)  in 
’,num2str(time_counter),'  sec’]), 
ylabel(’w_x(rad/sec)'),grid 

subplot(412),plot(time,w2),ylabel('w_y(rad/sec)'),ylabel('w_y(rad/sec)'),grid 

subplot(413),plot(time,w3),ylabel('w_z(rad/sec)'),ylabel('w_z(rad/sec)'),grid 

subplot(414),plot(time,w4),ylabel('w_t_o_t_a_l'),ylabel('w_t_o_t_a_l'),grid 

xlabel(’time(sec)') 

figure 

subplot(51  l),plot(time,pwr_u),title( [’Total  Power  =  ’,num2str(pwr3),’ 
W-hrs’]), ylabelC  W-hrs'),grid 

subplot(5 12), plot(time,watt_hrs_total), ylabelC  W-hrs'), ylabelC  W-hrs'), grid 
subplot(5 13),plot(time,ml  ,'r'),ylabel(’AmA2'),  ylabelC AmA2'), grid 
subplot(514),plot(time,m2,’g'),ylabel(’AmA2’),  ylabel(’AmA2’),grid 
subplot(5 1 5), plot(time, m3, 'b'), ylabel(' AmA2'),  ylabel(’AmA2’),grid 
xlabel('time(sec)') 

figure 

subplot(51  l),plot(time,pwrSS),title(['SS  Power  =  ’,num2str(pwr6),' 

W-hrs’]), ylabelC  W-hrs'), grid 

subplot(5 1 2), plot( time, pwr6), ylabelC  W-hrs'), ylabel(’  W-hrs'), grid 
subplot(5 13),plot(time,ml  ,’r'),ylabel(’AmA2'),  ylabel(’AmA2'),grid 
subplot(51 4), plot(time, m2, 'g'), ylabelC AmA2’),  ylabelC AmA2’), grid 
subplot(5 1 5), plot(time, m3, 'b'), ylabel(' AmA2'),  ylabelC AmA2’), grid 
xlabel(’time(sec)') 
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figure 

subplot(41  l),plot(time,rod3tot),title(['3  Rod  torque  =  ’,num2str(rod3tot),'  AmA2’]), 
ylabel('  AmA2'),grid 

subplot(4 12),plot(time,ml  ,'r'),ylabel(’AmA2'),  ylabel(’AmA2'),grid 
subplot(4 1 3),plot(time,m2,'g'),ylabel(’AmA2’),  ylabel(’AmA2'),grid 
subplot(414),plot(time,m3,’b’),ylabel(’AmA2’),  ylabel(’AmA2'),grid 
xlabel(’time(sec)’) 

figure 

plot(time,rod3tot),title(['3  Rod  torque  =  ’,num2str(rod3tot),’  AmA2’]), 

ylabel(’  AmA2'),grid 
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